From cfef013047a1c88788f8fc13eb33adf1000da8db Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Fri, 22 Mar 2019 10:09:47 +0100 Subject: [PATCH] Add tests for AD framework: ADOL-C --- tests/adolc/helper_retaping_01.cc | 297 +++ tests/adolc/helper_retaping_01.output | 3 + tests/adolc/helper_retaping_02.cc | 300 +++ tests/adolc/helper_retaping_02.output | 3 + tests/adolc/helper_retaping_03.cc | 304 +++ tests/adolc/helper_retaping_03.output | 35 + ...helper_scalar_coupled_2_components_01_1.cc | 39 + ...er_scalar_coupled_2_components_01_1.output | 3 + ...helper_scalar_coupled_2_components_01_2.cc | 39 + ...er_scalar_coupled_2_components_01_2.output | 3 + ...helper_scalar_coupled_2_components_02_1.cc | 39 + ...er_scalar_coupled_2_components_02_1.output | 3 + ...helper_scalar_coupled_2_components_02_2.cc | 39 + ...er_scalar_coupled_2_components_02_2.output | 3 + ...helper_scalar_coupled_2_components_03_1.cc | 39 + ...er_scalar_coupled_2_components_03_1.output | 3 + ...helper_scalar_coupled_2_components_03_2.cc | 39 + ...er_scalar_coupled_2_components_03_2.output | 3 + ...helper_scalar_coupled_2_components_04_1.cc | 39 + ...er_scalar_coupled_2_components_04_1.output | 3 + ...helper_scalar_coupled_2_components_04_2.cc | 39 + ...er_scalar_coupled_2_components_04_2.output | 3 + ...helper_scalar_coupled_2_components_05_1.cc | 39 + ...er_scalar_coupled_2_components_05_1.output | 3 + ...helper_scalar_coupled_2_components_05_2.cc | 39 + ...er_scalar_coupled_2_components_05_2.output | 3 + ...helper_scalar_coupled_2_components_06_1.cc | 43 + ...er_scalar_coupled_2_components_06_1.output | 3 + ...helper_scalar_coupled_2_components_06_2.cc | 43 + ...er_scalar_coupled_2_components_06_2.output | 3 + ...helper_scalar_coupled_2_components_07_1.cc | 43 + ...er_scalar_coupled_2_components_07_1.output | 3 + ...helper_scalar_coupled_2_components_07_2.cc | 43 + ...er_scalar_coupled_2_components_07_2.output | 3 + ...helper_scalar_coupled_2_components_08_1.cc | 45 + ...er_scalar_coupled_2_components_08_1.output | 3 + ...helper_scalar_coupled_2_components_08_2.cc | 45 + ...er_scalar_coupled_2_components_08_2.output | 3 + ...helper_scalar_coupled_3_components_01_1.cc | 44 + ...er_scalar_coupled_3_components_01_1.output | 3 + ...ponents_01_1.with_mpi=true.mpirun=2.output | 6 + ...helper_scalar_coupled_3_components_01_2.cc | 44 + ...er_scalar_coupled_3_components_01_2.output | 3 + ...ponents_01_2.with_mpi=true.mpirun=2.output | 6 + ...coupled_3_components_01_tapeless_only_2.cc | 48 + ...led_3_components_01_tapeless_only_2.output | 3 + tests/adolc/helper_scalar_failure_01.cc | 121 + tests/adolc/helper_scalar_failure_01.output | 2 + tests/adolc/helper_scalar_failure_02.cc | 123 + tests/adolc/helper_scalar_failure_02.output | 2 + tests/adolc/helper_scalar_failure_03.cc | 124 + tests/adolc/helper_scalar_failure_03.output | 2 + tests/adolc/helper_scalar_failure_04.cc | 123 + tests/adolc/helper_scalar_failure_04.output | 2 + tests/adolc/helper_scalar_failure_05.cc | 124 + tests/adolc/helper_scalar_failure_05.output | 2 + .../helper_scalar_single_component_01_1.cc | 38 + ...helper_scalar_single_component_01_1.output | 3 + .../helper_scalar_single_component_01_2.cc | 38 + ...helper_scalar_single_component_01_2.output | 3 + .../helper_scalar_single_component_02_1.cc | 38 + ...helper_scalar_single_component_02_1.output | 3 + .../helper_scalar_single_component_02_2.cc | 38 + ...helper_scalar_single_component_02_2.output | 3 + .../helper_scalar_single_component_03_1.cc | 38 + ...helper_scalar_single_component_03_1.output | 3 + .../helper_scalar_single_component_03_2.cc | 38 + ...helper_scalar_single_component_03_2.output | 3 + .../helper_scalar_single_component_04_1.cc | 39 + ...helper_scalar_single_component_04_1.output | 3 + .../helper_scalar_single_component_04_2.cc | 39 + ...helper_scalar_single_component_04_2.output | 3 + tests/adolc/helper_tape_index_01.cc | 338 +++ tests/adolc/helper_tape_index_01.output | 5 + .../helper_vector_1_indep_1_dep_vars_01_1.cc | 40 + ...lper_vector_1_indep_1_dep_vars_01_1.output | 3 + .../helper_vector_1_indep_1_dep_vars_01_2.cc | 40 + ...lper_vector_1_indep_1_dep_vars_01_2.output | 3 + .../helper_vector_1_indep_2_dep_vars_01_1.cc | 39 + ...lper_vector_1_indep_2_dep_vars_01_1.output | 3 + .../helper_vector_1_indep_2_dep_vars_01_2.cc | 39 + ...lper_vector_1_indep_2_dep_vars_01_2.output | 3 + .../helper_vector_2_indep_1_dep_vars_01_1.cc | 39 + ...lper_vector_2_indep_1_dep_vars_01_1.output | 3 + .../helper_vector_2_indep_1_dep_vars_01_2.cc | 39 + ...lper_vector_2_indep_1_dep_vars_01_2.output | 3 + .../helper_vector_2_indep_2_dep_vars_01_1.cc | 39 + ...lper_vector_2_indep_2_dep_vars_01_1.output | 3 + .../helper_vector_2_indep_2_dep_vars_01_2.cc | 39 + ...lper_vector_2_indep_2_dep_vars_01_2.output | 3 + tests/adolc/step-44-helper_scalar_01_1.cc | 68 + tests/adolc/step-44-helper_scalar_01_1.output | 11 + tests/adolc/step-44-helper_vector_01_1.cc | 68 + tests/adolc/step-44-helper_vector_01_1.output | 11 + tests/adolc/step-44-helper_vector_01_2.cc | 68 + tests/adolc/step-44-helper_vector_01_2.output | 11 + tests/adolc/step-44.cc | 2128 +++++++++++++++++ tests/adolc/step-44.output | 11 + .../adolc/symmetric_tensor_functions_02_1.cc | 40 + .../symmetric_tensor_functions_02_1.output | 3 + .../adolc/symmetric_tensor_functions_02_2.cc | 40 + .../symmetric_tensor_functions_02_2.output | 3 + .../adolc/symmetric_tensor_functions_03_1.cc | 40 + ...adolc_with_advanced_branching=false.output | 3 + ..._advanced_branching=true.expect=run.output | 1 + .../adolc/symmetric_tensor_functions_03_2.cc | 40 + .../symmetric_tensor_functions_03_2.output | 3 + .../adolc/symmetric_tensor_functions_04_1.cc | 40 + ...adolc_with_advanced_branching=false.output | 3 + ..._advanced_branching=true.expect=run.output | 1 + .../adolc/symmetric_tensor_functions_04_2.cc | 40 + .../symmetric_tensor_functions_04_2.output | 3 + tests/adolc/tensor_functions_02_1.cc | 40 + tests/adolc/tensor_functions_02_1.output | 3 + tests/adolc/tensor_functions_02_2.cc | 40 + tests/adolc/tensor_functions_02_2.output | 3 + tests/adolc/tensor_ops_01.cc | 236 ++ ...adolc_with_advanced_branching=false.output | 2 + ...ced_branching=true.expect=run.debug.output | 2 + ...ith_advanced_branching=true.release.output | 2 + tests/adolc/tensor_ops_02.cc | 105 + tests/adolc/tensor_ops_02.output | 2 + 122 files changed, 6506 insertions(+) create mode 100644 tests/adolc/helper_retaping_01.cc create mode 100644 tests/adolc/helper_retaping_01.output create mode 100644 tests/adolc/helper_retaping_02.cc create mode 100644 tests/adolc/helper_retaping_02.output create mode 100644 tests/adolc/helper_retaping_03.cc create mode 100644 tests/adolc/helper_retaping_03.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_01_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_01_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_01_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_01_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_02_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_02_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_02_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_02_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_03_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_03_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_03_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_03_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_04_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_04_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_04_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_04_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_05_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_05_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_05_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_05_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_06_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_06_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_06_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_06_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_07_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_07_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_07_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_07_2.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_08_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_08_1.output create mode 100644 tests/adolc/helper_scalar_coupled_2_components_08_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_2_components_08_2.output create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_1.cc create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_1.output create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_2.output create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc create mode 100644 tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output create mode 100644 tests/adolc/helper_scalar_failure_01.cc create mode 100644 tests/adolc/helper_scalar_failure_01.output create mode 100644 tests/adolc/helper_scalar_failure_02.cc create mode 100644 tests/adolc/helper_scalar_failure_02.output create mode 100644 tests/adolc/helper_scalar_failure_03.cc create mode 100644 tests/adolc/helper_scalar_failure_03.output create mode 100644 tests/adolc/helper_scalar_failure_04.cc create mode 100644 tests/adolc/helper_scalar_failure_04.output create mode 100644 tests/adolc/helper_scalar_failure_05.cc create mode 100644 tests/adolc/helper_scalar_failure_05.output create mode 100644 tests/adolc/helper_scalar_single_component_01_1.cc create mode 100644 tests/adolc/helper_scalar_single_component_01_1.output create mode 100644 tests/adolc/helper_scalar_single_component_01_2.cc create mode 100644 tests/adolc/helper_scalar_single_component_01_2.output create mode 100644 tests/adolc/helper_scalar_single_component_02_1.cc create mode 100644 tests/adolc/helper_scalar_single_component_02_1.output create mode 100644 tests/adolc/helper_scalar_single_component_02_2.cc create mode 100644 tests/adolc/helper_scalar_single_component_02_2.output create mode 100644 tests/adolc/helper_scalar_single_component_03_1.cc create mode 100644 tests/adolc/helper_scalar_single_component_03_1.output create mode 100644 tests/adolc/helper_scalar_single_component_03_2.cc create mode 100644 tests/adolc/helper_scalar_single_component_03_2.output create mode 100644 tests/adolc/helper_scalar_single_component_04_1.cc create mode 100644 tests/adolc/helper_scalar_single_component_04_1.output create mode 100644 tests/adolc/helper_scalar_single_component_04_2.cc create mode 100644 tests/adolc/helper_scalar_single_component_04_2.output create mode 100644 tests/adolc/helper_tape_index_01.cc create mode 100644 tests/adolc/helper_tape_index_01.output create mode 100644 tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc create mode 100644 tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output create mode 100644 tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc create mode 100644 tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output create mode 100644 tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc create mode 100644 tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output create mode 100644 tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc create mode 100644 tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output create mode 100644 tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc create mode 100644 tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output create mode 100644 tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc create mode 100644 tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output create mode 100644 tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc create mode 100644 tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output create mode 100644 tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc create mode 100644 tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output create mode 100644 tests/adolc/step-44-helper_scalar_01_1.cc create mode 100644 tests/adolc/step-44-helper_scalar_01_1.output create mode 100644 tests/adolc/step-44-helper_vector_01_1.cc create mode 100644 tests/adolc/step-44-helper_vector_01_1.output create mode 100644 tests/adolc/step-44-helper_vector_01_2.cc create mode 100644 tests/adolc/step-44-helper_vector_01_2.output create mode 100644 tests/adolc/step-44.cc create mode 100644 tests/adolc/step-44.output create mode 100644 tests/adolc/symmetric_tensor_functions_02_1.cc create mode 100644 tests/adolc/symmetric_tensor_functions_02_1.output create mode 100644 tests/adolc/symmetric_tensor_functions_02_2.cc create mode 100644 tests/adolc/symmetric_tensor_functions_02_2.output create mode 100644 tests/adolc/symmetric_tensor_functions_03_1.cc create mode 100644 tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output create mode 100644 tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output create mode 100644 tests/adolc/symmetric_tensor_functions_03_2.cc create mode 100644 tests/adolc/symmetric_tensor_functions_03_2.output create mode 100644 tests/adolc/symmetric_tensor_functions_04_1.cc create mode 100644 tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output create mode 100644 tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output create mode 100644 tests/adolc/symmetric_tensor_functions_04_2.cc create mode 100644 tests/adolc/symmetric_tensor_functions_04_2.output create mode 100644 tests/adolc/tensor_functions_02_1.cc create mode 100644 tests/adolc/tensor_functions_02_1.output create mode 100644 tests/adolc/tensor_functions_02_2.cc create mode 100644 tests/adolc/tensor_functions_02_2.output create mode 100644 tests/adolc/tensor_ops_01.cc create mode 100644 tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output create mode 100644 tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output create mode 100644 tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output create mode 100644 tests/adolc/tensor_ops_02.cc create mode 100644 tests/adolc/tensor_ops_02.output diff --git a/tests/adolc/helper_retaping_01.cc b/tests/adolc/helper_retaping_01.cc new file mode 100644 index 0000000000..a119926d5e --- /dev/null +++ b/tests/adolc/helper_retaping_01.cc @@ -0,0 +1,297 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check that a tape can be recorded and reused. In particular, we check +// that the state of the ADHelper can flipped from "record" to "reuse" +// and vice versa at any point in a simulation. Note that here the ADHelper +// will not have a persistent lifetime, so we expect that ADOL-C will keep +// any stored tape data as this object goes out of scope and is reconstructed +// at a later evaluation stage. +// +// Based off of helper_scalar_coupled_2_components_01.h, this test the +// evaluation of a coupled system (scalar + scalar components) +// using a helper class. +// +// AD number type: ADOL-C taped + +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +// Function and its derivatives +template +struct FunctionsTestScalarScalarCoupled +{ + static NumberType + psi(const NumberType &s1, const NumberType &s2) + { + return 2.0 * std::pow(s1, 4) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds1(const NumberType &s1, const NumberType &s2) + { + return 8.0 * std::pow(s1, 3) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds2(const NumberType &s1, const NumberType &s2) + { + return 6.0 * std::pow(s1, 4) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 2) * std::pow(s2, 3); + }; + + static NumberType + d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 3) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2) + { + return d2psi_ds2_ds1(s1, s2); + }; + + static NumberType + d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2) + { + return 12.0 * std::pow(s1, 4) * std::pow(s2, 1); + }; +}; + +template +void +test_scalar_scalar_coupled_one_run(const number_t s1, + const number_t s2, + const bool overwrite_tape) +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + typedef typename ADHelper::scalar_type ScalarNumberType; + + std::cout << "*** Test variables: Scalar + Scalar (coupled), " + << "dim = " << Utilities::to_string(dim) << ", " + << "Type code: " << static_cast(ad_type_code) << ". " + << "Overwrite tape: " << std::boolalpha << overwrite_tape + << std::endl; + + // Values computed from the AD energy function + ScalarNumberType psi; + Vector Dpsi; + FullMatrix D2psi; + + // Function and its derivatives + typedef FunctionsTestScalarScalarCoupled func_ad; + + // Setup the variable components and choose a value at which to + // evaluate the tape + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + const unsigned int n_AD_components = 2; + ADHelper ad_helper(n_AD_components); + if (overwrite_tape) + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + overwrite_tape, + true /*keep*/); + if (is_recording == true) + { + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi(func_ad::psi(s1_ad, s2_ad)); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Taped data..." << std::endl; + std::cout << "independent variable values: " << std::flush; + ad_helper.print_values(std::cout); + std::cout << "s1_ad: " << s1_ad << std::endl; + std::cout << "s2_ad: " << s2_ad << std::endl; + std::cout << "psi: " << psi << std::endl; + std::cout << std::endl; + } + else + { + std::cout + << "Using tape with different values for independent variables..." + << std::endl; + // Set a new evaluation point + ad_helper.activate_recorded_tape(tape_no); + ad_helper.set_independent_variable(s1, s1_dof); + ad_helper.set_independent_variable(s2, s2_dof); + } + + std::cout << "independent variable values: " << std::flush; + ad_helper.print_values(std::cout); + + // Compute the function value, gradient and hessian for the new evaluation + // point + psi = ad_helper.compute_value(); + ad_helper.compute_gradient(Dpsi); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + ad_helper.compute_hessian(D2psi); + } + + // Output the full stored function, gradient vector and hessian matrix + std::cout << "psi: " << psi << std::endl; + std::cout << "Dpsi: \n"; + Dpsi.print(std::cout); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + std::cout << "D2psi: \n"; + D2psi.print_formatted(std::cout, 3, true, 0, "0.0"); + } + + // Extract components of the solution + const ScalarNumberType dpsi_ds1 = + ad_helper.extract_gradient_component(Dpsi, s1_dof); + const ScalarNumberType dpsi_ds2 = + ad_helper.extract_gradient_component(Dpsi, s2_dof); + std::cout << "extracted Dpsi (s1): " << dpsi_ds1 << "\n" + << "extracted Dpsi (s2): " << dpsi_ds2 << "\n"; + + // Verify the result + typedef FunctionsTestScalarScalarCoupled func; + static const ScalarNumberType tol = + 1e8 * std::numeric_limits::epsilon(); + std::cout << "dpsi_ds1: " << dpsi_ds1 << std::endl; + std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1, s2) << std::endl; + std::cout << "diff: " << std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) + << std::endl; + std::cout << "dpsi_ds2: " << dpsi_ds2 << std::endl; + std::cout << "func::dpsi_ds2(s1,s2): " << func::dpsi_ds2(s1, s2) << std::endl; + std::cout << "diff: " << std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) + << std::endl; + Assert(std::abs(psi - func::psi(s1, s2)) < tol, + ExcMessage("No match for function value.")); + Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + const ScalarNumberType d2psi_ds1_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof); + const ScalarNumberType d2psi_ds1_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof); + std::cout << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n" + << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n" + << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n" + << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n" + << std::endl; + std::cout << "d2psi_ds1_ds1: " << d2psi_ds1_ds1 << std::endl; + std::cout << "func::d2psi_ds1_ds1(s1,s2): " << func::d2psi_ds1_ds1(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) + << std::endl; + std::cout << "d2psi_ds2_ds1: " << d2psi_ds2_ds1 << std::endl; + std::cout << "func::d2psi_ds2_ds1(s1,s2): " << func::d2psi_ds2_ds1(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) + << std::endl; + std::cout << "d2psi_ds1_ds2: " << d2psi_ds1_ds2 << std::endl; + std::cout << "func::d2psi_ds1_ds2(s1,s2): " << func::d2psi_ds1_ds2(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) + << std::endl; + std::cout << "d2psi_ds2_ds2: " << d2psi_ds2_ds2 << std::endl; + std::cout << "func::d2psi_ds2_ds2(s1,s2): " << func::d2psi_ds2_ds2(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) + << std::endl; + Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + } +} + + + +template +void +test_scalar_scalar_coupled(const unsigned int n_runs, + const unsigned int retaping_modulo) +{ + for (unsigned int r = 0; r < n_runs; ++r) + { + const bool overwrite_tape = (r % retaping_modulo == 0); + const number_t s1 = 3.1 * (1.0 + r / 10.0); + const number_t s2 = 5.9 * (1.0 + r / 10.0); + test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>( + s1, s2, overwrite_tape); + } +} + + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(12, 4); + test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(4, 2); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_retaping_01.output b/tests/adolc/helper_retaping_01.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_retaping_01.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_retaping_02.cc b/tests/adolc/helper_retaping_02.cc new file mode 100644 index 0000000000..485540160b --- /dev/null +++ b/tests/adolc/helper_retaping_02.cc @@ -0,0 +1,300 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check that a tape can be recorded and reused. In particular, we check +// that the state of the ADHelper can flipped from "record" to "reuse" +// and vice versa at any point in a simulation. As opposed to +// helper_retaping_01, the ADHelper is now a static object (i.e. it has +// persistent lifetime, as might be the case if an ADHelper is stored +// as quadrature point member data). +// +// Based off of helper_scalar_coupled_2_components_01.h, this test the +// evaluation of a coupled system (scalar + scalar components) +// using a helper class. +// +// AD number type: ADOL-C taped + +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +// Function and its derivatives +template +struct FunctionsTestScalarScalarCoupled +{ + static NumberType + psi(const NumberType &s1, const NumberType &s2) + { + return 2.0 * std::pow(s1, 4) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds1(const NumberType &s1, const NumberType &s2) + { + return 8.0 * std::pow(s1, 3) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds2(const NumberType &s1, const NumberType &s2) + { + return 6.0 * std::pow(s1, 4) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 2) * std::pow(s2, 3); + }; + + static NumberType + d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 3) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2) + { + return d2psi_ds2_ds1(s1, s2); + }; + + static NumberType + d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2) + { + return 12.0 * std::pow(s1, 4) * std::pow(s2, 1); + }; +}; + +template +void +test_scalar_scalar_coupled_one_run(const number_t s1, + const number_t s2, + const bool overwrite_tape) +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + typedef typename ADHelper::scalar_type ScalarNumberType; + + std::cout << "*** Test variables: Scalar + Scalar (coupled), " + << "dim = " << Utilities::to_string(dim) << ", " + << "Type code: " << static_cast(ad_type_code) << ". " + << "Overwrite tape: " << std::boolalpha << overwrite_tape + << std::endl; + + // Values computed from the AD energy function + ScalarNumberType psi; + Vector Dpsi; + FullMatrix D2psi; + + // Function and its derivatives + typedef FunctionsTestScalarScalarCoupled func_ad; + + // Setup the variable components and choose a value at which to + // evaluate the tape + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + const unsigned int n_AD_components = 2; + static ADHelper ad_helper(n_AD_components); + if (overwrite_tape) + { + ad_helper.reset(); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + } + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + overwrite_tape, + true /*keep*/); + if (is_recording == true) + { + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi(func_ad::psi(s1_ad, s2_ad)); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Taped data..." << std::endl; + std::cout << "independent variable values: " << std::flush; + ad_helper.print_values(std::cout); + std::cout << "s1_ad: " << s1_ad << std::endl; + std::cout << "s2_ad: " << s2_ad << std::endl; + std::cout << "psi: " << psi << std::endl; + std::cout << std::endl; + } + else + { + std::cout + << "Using tape with different values for independent variables..." + << std::endl; + // Set a new evaluation point + ad_helper.activate_recorded_tape(tape_no); + ad_helper.set_independent_variable(s1, s1_dof); + ad_helper.set_independent_variable(s2, s2_dof); + } + + std::cout << "independent variable values: " << std::flush; + ad_helper.print_values(std::cout); + + // Compute the function value, gradient and hessian for the new evaluation + // point + psi = ad_helper.compute_value(); + ad_helper.compute_gradient(Dpsi); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + ad_helper.compute_hessian(D2psi); + } + + // Output the full stored function, gradient vector and hessian matrix + std::cout << "psi: " << psi << std::endl; + std::cout << "Dpsi: \n"; + Dpsi.print(std::cout); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + std::cout << "D2psi: \n"; + D2psi.print_formatted(std::cout, 3, true, 0, "0.0"); + } + + // Extract components of the solution + const ScalarNumberType dpsi_ds1 = + ad_helper.extract_gradient_component(Dpsi, s1_dof); + const ScalarNumberType dpsi_ds2 = + ad_helper.extract_gradient_component(Dpsi, s2_dof); + std::cout << "extracted Dpsi (s1): " << dpsi_ds1 << "\n" + << "extracted Dpsi (s2): " << dpsi_ds2 << "\n"; + + // Verify the result + typedef FunctionsTestScalarScalarCoupled func; + static const ScalarNumberType tol = + 1e8 * std::numeric_limits::epsilon(); + std::cout << "dpsi_ds1: " << dpsi_ds1 << std::endl; + std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1, s2) << std::endl; + std::cout << "diff: " << std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) + << std::endl; + std::cout << "dpsi_ds2: " << dpsi_ds2 << std::endl; + std::cout << "func::dpsi_ds2(s1,s2): " << func::dpsi_ds2(s1, s2) << std::endl; + std::cout << "diff: " << std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) + << std::endl; + Assert(std::abs(psi - func::psi(s1, s2)) < tol, + ExcMessage("No match for function value.")); + Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + const ScalarNumberType d2psi_ds1_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof); + const ScalarNumberType d2psi_ds1_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof); + std::cout << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n" + << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n" + << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n" + << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n" + << std::endl; + std::cout << "d2psi_ds1_ds1: " << d2psi_ds1_ds1 << std::endl; + std::cout << "func::d2psi_ds1_ds1(s1,s2): " << func::d2psi_ds1_ds1(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) + << std::endl; + std::cout << "d2psi_ds2_ds1: " << d2psi_ds2_ds1 << std::endl; + std::cout << "func::d2psi_ds2_ds1(s1,s2): " << func::d2psi_ds2_ds1(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) + << std::endl; + std::cout << "d2psi_ds1_ds2: " << d2psi_ds1_ds2 << std::endl; + std::cout << "func::d2psi_ds1_ds2(s1,s2): " << func::d2psi_ds1_ds2(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) + << std::endl; + std::cout << "d2psi_ds2_ds2: " << d2psi_ds2_ds2 << std::endl; + std::cout << "func::d2psi_ds2_ds2(s1,s2): " << func::d2psi_ds2_ds2(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) + << std::endl; + Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + } +} + + + +template +void +test_scalar_scalar_coupled(const unsigned int n_runs, + const unsigned int retaping_modulo) +{ + for (unsigned int r = 0; r < n_runs; ++r) + { + const bool overwrite_tape = (r % retaping_modulo == 0); + const number_t s1 = 3.1 * (1.0 + r / 10.0); + const number_t s2 = 5.9 * (1.0 + r / 10.0); + test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>( + s1, s2, overwrite_tape); + } +} + + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(12, 4); + test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(4, 2); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_retaping_02.output b/tests/adolc/helper_retaping_02.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_retaping_02.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_retaping_03.cc b/tests/adolc/helper_retaping_03.cc new file mode 100644 index 0000000000..9a982257c8 --- /dev/null +++ b/tests/adolc/helper_retaping_03.cc @@ -0,0 +1,304 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check that a tape can be recorded and reused. Here we specifically check +// the status of the is_registered_tape() flag, and what happens after all +// tapes are cleared. We explicitly do not allow tapes to be overwritten. +// Note that this also checks if we can safely query the existence of a +// tape that has not yet been traced (i.e. it is non-existent). +// +// Based off of helper_scalar_coupled_2_components_01.h, this test the +// evaluation of a coupled system (scalar + scalar components) +// using a helper class. +// +// AD number type: ADOL-C taped + +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +// Function and its derivatives +template +struct FunctionsTestScalarScalarCoupled +{ + static NumberType + psi(const NumberType &s1, const NumberType &s2) + { + return 2.0 * std::pow(s1, 4) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds1(const NumberType &s1, const NumberType &s2) + { + return 8.0 * std::pow(s1, 3) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds2(const NumberType &s1, const NumberType &s2) + { + return 6.0 * std::pow(s1, 4) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 2) * std::pow(s2, 3); + }; + + static NumberType + d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 3) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2) + { + return d2psi_ds2_ds1(s1, s2); + }; + + static NumberType + d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2) + { + return 12.0 * std::pow(s1, 4) * std::pow(s2, 1); + }; +}; + +template +void +test_scalar_scalar_coupled_one_run(const number_t s1, + const number_t s2, + const bool clear_tapes) +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + typedef typename ADHelper::scalar_type ScalarNumberType; + + std::cout << "*** Test variables: Scalar + Scalar (coupled), " + << "dim = " << Utilities::to_string(dim) << ", " + << "Type code: " << static_cast(ad_type_code) << ". " + << "Clear tapes: " << std::boolalpha << clear_tapes << std::endl; + + // Values computed from the AD energy function + ScalarNumberType psi; + Vector Dpsi; + FullMatrix D2psi; + + // Function and its derivatives + typedef FunctionsTestScalarScalarCoupled func_ad; + + // Setup the variable components and choose a value at which to + // evaluate the tape + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + const unsigned int n_AD_components = 2; + ADHelper ad_helper(n_AD_components); + if (clear_tapes) + { + ad_helper.reset(); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + } + + // Configure tape + const int tape_no = 1; + deallog << "is_registered_tape? " << std::boolalpha + << ad_helper.is_registered_tape(tape_no) << std::endl; + + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + false /*overwrite_tape*/, + true /*keep*/); + if (is_recording == true) + { + deallog << "Recording new tape" << std::endl; + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi(func_ad::psi(s1_ad, s2_ad)); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Taped data..." << std::endl; + std::cout << "independent variable values: " << std::flush; + ad_helper.print_values(std::cout); + std::cout << "s1_ad: " << s1_ad << std::endl; + std::cout << "s2_ad: " << s2_ad << std::endl; + std::cout << "psi: " << psi << std::endl; + std::cout << std::endl; + } + else + { + deallog << "Reusing prerecorded tape" << std::endl; + std::cout + << "Using tape with different values for independent variables..." + << std::endl; + // Set a new evaluation point + ad_helper.activate_recorded_tape(tape_no); + ad_helper.set_independent_variable(s1, s1_dof); + ad_helper.set_independent_variable(s2, s2_dof); + } + + std::cout << "independent variable values: " << std::flush; + ad_helper.print_values(std::cout); + + // Compute the function value, gradient and hessian for the new evaluation + // point + psi = ad_helper.compute_value(); + ad_helper.compute_gradient(Dpsi); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + ad_helper.compute_hessian(D2psi); + } + + // Output the full stored function, gradient vector and hessian matrix + std::cout << "psi: " << psi << std::endl; + std::cout << "Dpsi: \n"; + Dpsi.print(std::cout); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + std::cout << "D2psi: \n"; + D2psi.print_formatted(std::cout, 3, true, 0, "0.0"); + } + + // Extract components of the solution + const ScalarNumberType dpsi_ds1 = + ad_helper.extract_gradient_component(Dpsi, s1_dof); + const ScalarNumberType dpsi_ds2 = + ad_helper.extract_gradient_component(Dpsi, s2_dof); + std::cout << "extracted Dpsi (s1): " << dpsi_ds1 << "\n" + << "extracted Dpsi (s2): " << dpsi_ds2 << "\n"; + + // Verify the result + typedef FunctionsTestScalarScalarCoupled func; + static const ScalarNumberType tol = + 1e8 * std::numeric_limits::epsilon(); + std::cout << "dpsi_ds1: " << dpsi_ds1 << std::endl; + std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1, s2) << std::endl; + std::cout << "diff: " << std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) + << std::endl; + std::cout << "dpsi_ds2: " << dpsi_ds2 << std::endl; + std::cout << "func::dpsi_ds2(s1,s2): " << func::dpsi_ds2(s1, s2) << std::endl; + std::cout << "diff: " << std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) + << std::endl; + Assert(std::abs(psi - func::psi(s1, s2)) < tol, + ExcMessage("No match for function value.")); + Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + const ScalarNumberType d2psi_ds1_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof); + const ScalarNumberType d2psi_ds1_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof); + std::cout << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n" + << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n" + << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n" + << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n" + << std::endl; + std::cout << "d2psi_ds1_ds1: " << d2psi_ds1_ds1 << std::endl; + std::cout << "func::d2psi_ds1_ds1(s1,s2): " << func::d2psi_ds1_ds1(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) + << std::endl; + std::cout << "d2psi_ds2_ds1: " << d2psi_ds2_ds1 << std::endl; + std::cout << "func::d2psi_ds2_ds1(s1,s2): " << func::d2psi_ds2_ds1(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) + << std::endl; + std::cout << "d2psi_ds1_ds2: " << d2psi_ds1_ds2 << std::endl; + std::cout << "func::d2psi_ds1_ds2(s1,s2): " << func::d2psi_ds1_ds2(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) + << std::endl; + std::cout << "d2psi_ds2_ds2: " << d2psi_ds2_ds2 << std::endl; + std::cout << "func::d2psi_ds2_ds2(s1,s2): " << func::d2psi_ds2_ds2(s1, s2) + << std::endl; + std::cout << "diff: " + << std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) + << std::endl; + Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + } +} + + + +template +void +test_scalar_scalar_coupled(const unsigned int n_runs, + const unsigned int clearing_modulo) +{ + for (unsigned int r = 0; r < n_runs; ++r) + { + const bool clear_tape = (r % clearing_modulo == 0); + const number_t s1 = 3.1 * (1.0 + r / 10.0); + const number_t s2 = 5.9 * (1.0 + r / 10.0); + test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>(s1, + s2, + clear_tape); + } +} + + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(12, 4); + test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(4, 2); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_retaping_03.output b/tests/adolc/helper_retaping_03.output new file mode 100644 index 0000000000..1dffec25c4 --- /dev/null +++ b/tests/adolc/helper_retaping_03.output @@ -0,0 +1,35 @@ + +DEAL:Double::is_registered_tape? false +DEAL:Double::Recording new tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? false +DEAL:Double::Recording new tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? false +DEAL:Double::Recording new tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? false +DEAL:Double::Recording new tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::is_registered_tape? false +DEAL:Double::Recording new tape +DEAL:Double::is_registered_tape? true +DEAL:Double::Reusing prerecorded tape +DEAL:Double::OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_1.cc b/tests/adolc/helper_scalar_coupled_2_components_01_1.cc new file mode 100644 index 0000000000..08cb486b9e --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_01_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (scalar + scalar components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(); + test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_1.output b/tests/adolc/helper_scalar_coupled_2_components_01_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_01_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_2.cc b/tests/adolc/helper_scalar_coupled_2_components_01_2.cc new file mode 100644 index 0000000000..e86a3c2031 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_01_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (scalar + scalar components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_tapeless>(); + test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_2.output b/tests/adolc/helper_scalar_coupled_2_components_01_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_01_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_1.cc b/tests/adolc/helper_scalar_coupled_2_components_02_1.cc new file mode 100644 index 0000000000..a7fc5036dd --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_02_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (vector + scalar components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_vector_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(); + test_vector_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_1.output b/tests/adolc/helper_scalar_coupled_2_components_02_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_02_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_2.cc b/tests/adolc/helper_scalar_coupled_2_components_02_2.cc new file mode 100644 index 0000000000..a5ddb71e45 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_02_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (vector + scalar components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_vector_scalar_coupled<2, double, AD::NumberTypes::adolc_tapeless>(); + test_vector_scalar_coupled<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_2.output b/tests/adolc/helper_scalar_coupled_2_components_02_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_02_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_1.cc b/tests/adolc/helper_scalar_coupled_2_components_03_1.cc new file mode 100644 index 0000000000..905a93f676 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_03_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + scalar components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_03.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(); + test_tensor_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_1.output b/tests/adolc/helper_scalar_coupled_2_components_03_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_03_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_2.cc b/tests/adolc/helper_scalar_coupled_2_components_03_2.cc new file mode 100644 index 0000000000..67b3d5a2c3 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_03_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + scalar components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_03.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor_scalar_coupled<2, double, AD::NumberTypes::adolc_tapeless>(); + test_tensor_scalar_coupled<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_2.output b/tests/adolc/helper_scalar_coupled_2_components_03_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_03_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_1.cc b/tests/adolc/helper_scalar_coupled_2_components_04_1.cc new file mode 100644 index 0000000000..f16578ae29 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_04_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (vector + vector components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_04.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_vector_vector_coupled<2, double, AD::NumberTypes::adolc_taped>(); + test_vector_vector_coupled<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_1.output b/tests/adolc/helper_scalar_coupled_2_components_04_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_04_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_2.cc b/tests/adolc/helper_scalar_coupled_2_components_04_2.cc new file mode 100644 index 0000000000..295128d357 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_04_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (vector + vector components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_04.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_vector_vector_coupled<2, double, AD::NumberTypes::adolc_tapeless>(); + test_vector_vector_coupled<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_2.output b/tests/adolc/helper_scalar_coupled_2_components_04_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_04_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_1.cc b/tests/adolc/helper_scalar_coupled_2_components_05_1.cc new file mode 100644 index 0000000000..bb60dc5429 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_05_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + vector components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_05.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor_vector_coupled<2, double, AD::NumberTypes::adolc_taped>(); + test_tensor_vector_coupled<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_1.output b/tests/adolc/helper_scalar_coupled_2_components_05_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_05_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_2.cc b/tests/adolc/helper_scalar_coupled_2_components_05_2.cc new file mode 100644 index 0000000000..eb197e8891 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_05_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + vector components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_05.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor_vector_coupled<2, double, AD::NumberTypes::adolc_tapeless>(); + test_tensor_vector_coupled<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_2.output b/tests/adolc/helper_scalar_coupled_2_components_05_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_05_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_1.cc b/tests/adolc/helper_scalar_coupled_2_components_06_1.cc new file mode 100644 index 0000000000..1c49818ce2 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_06_1.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (symmetric tensor + vector components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_06.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor_vector_coupled<2, + double, + AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor_vector_coupled<3, + double, + AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_1.output b/tests/adolc/helper_scalar_coupled_2_components_06_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_06_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_2.cc b/tests/adolc/helper_scalar_coupled_2_components_06_2.cc new file mode 100644 index 0000000000..9252106fd3 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_06_2.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (symmetric tensor + vector components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_06.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor_vector_coupled<2, + double, + AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor_vector_coupled<3, + double, + AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_2.output b/tests/adolc/helper_scalar_coupled_2_components_06_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_06_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_1.cc b/tests/adolc/helper_scalar_coupled_2_components_07_1.cc new file mode 100644 index 0000000000..c3aabfbf63 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_07_1.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (symmetric tensor + tensor components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_07.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor_tensor_coupled<2, + double, + AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor_tensor_coupled<3, + double, + AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_1.output b/tests/adolc/helper_scalar_coupled_2_components_07_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_07_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_2.cc b/tests/adolc/helper_scalar_coupled_2_components_07_2.cc new file mode 100644 index 0000000000..3b9dc6828b --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_07_2.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (symmetric tensor + tensor components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_07.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor_tensor_coupled<2, + double, + AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor_tensor_coupled<3, + double, + AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_2.output b/tests/adolc/helper_scalar_coupled_2_components_07_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_07_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_1.cc b/tests/adolc/helper_scalar_coupled_2_components_08_1.cc new file mode 100644 index 0000000000..126f69b731 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_08_1.cc @@ -0,0 +1,45 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (symmetric tensor + symmetric tensor +// components) using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_2_components_08.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor_symmetric_tensor_coupled< + 2, + double, + AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor_symmetric_tensor_coupled< + 3, + double, + AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_1.output b/tests/adolc/helper_scalar_coupled_2_components_08_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_08_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_2.cc b/tests/adolc/helper_scalar_coupled_2_components_08_2.cc new file mode 100644 index 0000000000..2bc24494ef --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_08_2.cc @@ -0,0 +1,45 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (symmetric tensor + symmetric tensor +// components) using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_2_components_08.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor_symmetric_tensor_coupled< + 2, + double, + AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor_symmetric_tensor_coupled< + 3, + double, + AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_2.output b/tests/adolc/helper_scalar_coupled_2_components_08_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_2_components_08_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_1.cc b/tests/adolc/helper_scalar_coupled_3_components_01_1.cc new file mode 100644 index 0000000000..db7bb573ab --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_1.cc @@ -0,0 +1,44 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + vector + scalar components) +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_coupled_3_components_01.h" +#include "../tests.h" + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll mpi_log; + + deallog.push("Double"); + { + test_tensor_vector_scalar_coupled<2, + double, + AD::NumberTypes::adolc_taped>(); + test_tensor_vector_scalar_coupled<3, + double, + AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_1.output b/tests/adolc/helper_scalar_coupled_3_components_01_1.output new file mode 100644 index 0000000000..aeccfbdd06 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_1.output @@ -0,0 +1,3 @@ + +DEAL:0:Double:OK +DEAL:0:OK diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output b/tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output new file mode 100644 index 0000000000..44481250fc --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output @@ -0,0 +1,6 @@ + +DEAL:0:Double::OK +DEAL:0:OK + +DEAL:1:Double::OK +DEAL:1:OK \ No newline at end of file diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_2.cc b/tests/adolc/helper_scalar_coupled_3_components_01_2.cc new file mode 100644 index 0000000000..8a5c950161 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_2.cc @@ -0,0 +1,44 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + vector + scalar components) +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_3_components_01.h" +#include "../tests.h" + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll mpi_log; + + deallog.push("Double"); + { + test_tensor_vector_scalar_coupled<2, + double, + AD::NumberTypes::adolc_tapeless>(); + test_tensor_vector_scalar_coupled<3, + double, + AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_2.output b/tests/adolc/helper_scalar_coupled_3_components_01_2.output new file mode 100644 index 0000000000..aeccfbdd06 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_2.output @@ -0,0 +1,3 @@ + +DEAL:0:Double:OK +DEAL:0:OK diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output b/tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output new file mode 100644 index 0000000000..44481250fc --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output @@ -0,0 +1,6 @@ + +DEAL:0:Double::OK +DEAL:0:OK + +DEAL:1:Double::OK +DEAL:1:OK \ No newline at end of file diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc b/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc new file mode 100644 index 0000000000..f13cecaee2 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc @@ -0,0 +1,48 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a coupled system (tensor + vector + scalar components) +// using a helper class. +// This test is based off of helper_scalar_coupled_3_components_01.h, and checks +// that everything still works in tapeless only mode (i.e. when the +// start_recording_operations and stop_recording_operations calls are +// removed). +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_coupled_3_components_01_tapeless_only.h" +#include "../tests.h" + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll mpi_log; + + deallog.push("Double"); + { + test_tensor_vector_scalar_coupled<2, + double, + AD::NumberTypes::adolc_tapeless>(); + test_tensor_vector_scalar_coupled<3, + double, + AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output b/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output new file mode 100644 index 0000000000..aeccfbdd06 --- /dev/null +++ b/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output @@ -0,0 +1,3 @@ + +DEAL:0:Double:OK +DEAL:0:OK diff --git a/tests/adolc/helper_scalar_failure_01.cc b/tests/adolc/helper_scalar_failure_01.cc new file mode 100644 index 0000000000..09aa0b72c5 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_01.cc @@ -0,0 +1,121 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test that some cases that we expect to fail do in fact do so: +// Overlapping indices of independent variables + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +template +bool +test_overlapping_ind_var_components() +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + + std::cout << "*** Test variables: Scalar + Scalar (coupled), " + << (AD::ADNumberTraits::is_taped == true ? "Taped" : + "Tapeless") + << std::endl; + + try + { + // NOTE: Both extractors set to same index + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(0); + const unsigned int n_AD_components = 2; + ADHelper ad_helper(n_AD_components); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + + double s1 = 3.1; + double s2 = 5.9; + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + true /*overwrite_tape*/, + true /*keep*/); + + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi(s1_ad * s2_ad); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Unexpected success!" << std::endl; + return true; + } + catch (...) + { + std::cout << "Expected failure." << std::endl; + return false; + } +} + +int +main() +{ + initlog(); + + deal_II_exceptions::disable_abort_on_exception(); +#if DEBUG + // Asserts should be triggered + const bool expected_result = false; +#else + // User beware: Asserts ignored + const bool expected_result = true; +#endif + + const unsigned int dim = 2; + AssertThrow( + (test_overlapping_ind_var_components() == + expected_result), + ExcInternalError()); + AssertThrow( + (test_overlapping_ind_var_components() == + expected_result), + ExcInternalError()); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_failure_01.output b/tests/adolc/helper_scalar_failure_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_01.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/helper_scalar_failure_02.cc b/tests/adolc/helper_scalar_failure_02.cc new file mode 100644 index 0000000000..59757f1d27 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_02.cc @@ -0,0 +1,123 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test that some cases that we expect to fail do in fact do so: +// Overlapping indices of independent variables + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +template +bool +test_overlapping_ind_var_components() +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + + std::cout << "*** Test variables: Vector + Scalar (coupled), " + << (AD::ADNumberTraits::is_taped == true ? "Taped" : + "Tapeless") + << std::endl; + + try + { + // NOTE: Second extractor overlaps first ont. + const FEValuesExtractors::Vector v1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + const unsigned int n_AD_components = + Tensor<1, dim>::n_independent_components + 1; + ADHelper ad_helper(n_AD_components); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + + Tensor<1, dim> v1; + double s2 = 5.9; + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + true /*overwrite_tape*/, + true /*keep*/); + + ad_helper.register_independent_variable(v1, v1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const Tensor<1, dim, ADNumberType> v1_ad = + ad_helper.get_sensitive_variables(v1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi((v1_ad * v1_ad) * s2_ad); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Unexpected success!" << std::endl; + return true; + } + catch (...) + { + std::cout << "Expected failure." << std::endl; + return false; + } +} + +int +main() +{ + initlog(); + + deal_II_exceptions::disable_abort_on_exception(); +#if DEBUG + // Asserts should be triggered + const bool expected_result = false; +#else + // User beware: Asserts ignored + const bool expected_result = true; +#endif + + const unsigned int dim = 2; + AssertThrow( + (test_overlapping_ind_var_components() == + expected_result), + ExcInternalError()); + AssertThrow( + (test_overlapping_ind_var_components() == + expected_result), + ExcInternalError()); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_failure_02.output b/tests/adolc/helper_scalar_failure_02.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_02.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/helper_scalar_failure_03.cc b/tests/adolc/helper_scalar_failure_03.cc new file mode 100644 index 0000000000..110c4e04e6 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_03.cc @@ -0,0 +1,124 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test that some cases that we expect to fail do in fact do so: +// Overlapping indices of independent variables + +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +template +bool +test_overlapping_ind_var_components() +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + + std::cout << "*** Test variables: Tensor + Vector (coupled), " + << (AD::ADNumberTraits::is_taped == true ? "Taped" : + "Tapeless") + << std::endl; + + try + { + // NOTE: Second extractor overlaps first ont. + const FEValuesExtractors::Tensor<2> t1_dof(0); + const FEValuesExtractors::Vector v2_dof(dim + 1); + const unsigned int n_AD_components = + Tensor<2, dim>::n_independent_components + + Tensor<1, dim>::n_independent_components; + ADHelper ad_helper(n_AD_components); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + + Tensor<2, dim> t1; + Tensor<1, dim> v2; + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + true /*overwrite_tape*/, + true /*keep*/); + + ad_helper.register_independent_variable(t1, t1_dof); + ad_helper.register_independent_variable(v2, v2_dof); + + const Tensor<2, dim, ADNumberType> t1_ad = + ad_helper.get_sensitive_variables(t1_dof); + const Tensor<1, dim, ADNumberType> v2_ad = + ad_helper.get_sensitive_variables(v2_dof); + + const ADNumberType psi((v2_ad * v2_ad) * trace(t1_ad)); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Unexpected success!" << std::endl; + return true; + } + catch (...) + { + std::cout << "Expected failure." << std::endl; + return false; + } +} + +int +main() +{ + initlog(); + + deal_II_exceptions::disable_abort_on_exception(); +#if DEBUG + // Asserts should be triggered + const bool expected_result = false; +#else + // User beware: Asserts ignored + const bool expected_result = true; +#endif + + const unsigned int dim = 2; + AssertThrow( + (test_overlapping_ind_var_components() == + expected_result), + ExcInternalError()); + AssertThrow( + (test_overlapping_ind_var_components() == + expected_result), + ExcInternalError()); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_failure_03.output b/tests/adolc/helper_scalar_failure_03.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_03.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/helper_scalar_failure_04.cc b/tests/adolc/helper_scalar_failure_04.cc new file mode 100644 index 0000000000..bd757b8991 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_04.cc @@ -0,0 +1,123 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test that some cases that we expect to fail do in fact do so: +// Incorrect setting of number of independent variables + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +template +bool +test_incorrect_no_of_ind_vars() +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + + std::cout << "*** Test variables: Scalar + Scalar (coupled), " + << (AD::ADNumberTraits::is_taped == true ? "Taped" : + "Tapeless") + << std::endl; + + try + { + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + // NOTE: Only 1 independent variable set instead of 2 + const unsigned int n_AD_components = 1; + ADHelper ad_helper(n_AD_components); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + + double s1 = 3.1; + double s2 = 5.9; + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + true /*overwrite_tape*/, + true /*keep*/); + + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables( + s2_dof); // This line might lead to memory corruption in tapeless mode. + + const ADNumberType psi(s1_ad * s2_ad); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Unexpected success!" << std::endl; + return true; + } + catch (...) + { + std::cout << "Expected failure." << std::endl; + return false; + } +} + +int +main() +{ + initlog(); + + deal_II_exceptions::disable_abort_on_exception(); +#if DEBUG + // Asserts should be triggered + const bool expected_result_taped = false; + const bool expected_result_tapeless = false; +#else + // User beware: Asserts ignored + const bool expected_result_taped = true; + const bool expected_result_tapeless = false; +#endif + + const unsigned int dim = 2; + AssertThrow((test_incorrect_no_of_ind_vars() == + expected_result_taped), + ExcInternalError()); + AssertThrow( + (test_incorrect_no_of_ind_vars() == + expected_result_tapeless), + ExcInternalError()); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_failure_04.output b/tests/adolc/helper_scalar_failure_04.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_04.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/helper_scalar_failure_05.cc b/tests/adolc/helper_scalar_failure_05.cc new file mode 100644 index 0000000000..63e41a6a80 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_05.cc @@ -0,0 +1,124 @@ + +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test that some cases that we expect to fail do in fact do so: +// Registering more than one dependent variable + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +template +bool +test_register_multiple_dep_vars() +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + + std::cout << "*** Test variables: Scalar + Scalar (coupled), " + << (AD::ADNumberTraits::is_taped == true ? "Taped" : + "Tapeless") + << std::endl; + + try + { + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + const unsigned int n_AD_components = 2; + ADHelper ad_helper(n_AD_components); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + + double s1 = 3.1; + double s2 = 5.9; + + // Configure tape + const int tape_no = 1; + const bool is_recording = + ad_helper.start_recording_operations(tape_no /*material_id*/, + true /*overwrite_tape*/, + true /*keep*/); + + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi(s1_ad * s2_ad); + + // NOTE: Should only be calling this function once for this helper class + ad_helper.register_dependent_variable(psi); + ad_helper.register_dependent_variable(psi); // Wrong! + + ad_helper.stop_recording_operations(false /*write_tapes_to_file*/); + + std::cout << "Unexpected success!" << std::endl; + return true; + } + catch (...) + { + std::cout << "Expected failure." << std::endl; + return false; + } +} + +int +main() +{ + initlog(); + + deal_II_exceptions::disable_abort_on_exception(); +#if DEBUG + // Asserts should be triggered + const bool expected_result = false; +#else + // User beware: Asserts ignored + const bool expected_result = true; +#endif + + const unsigned int dim = 2; + AssertThrow( + (test_register_multiple_dep_vars() == + expected_result), + ExcInternalError()); + AssertThrow( + (test_register_multiple_dep_vars() == + expected_result), + ExcInternalError()); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_failure_05.output b/tests/adolc/helper_scalar_failure_05.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/helper_scalar_failure_05.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_01_1.cc b/tests/adolc/helper_scalar_single_component_01_1.cc new file mode 100644 index 0000000000..0ef35c8901 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_01_1.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (scalar) system using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_single_component_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar<2, double, AD::NumberTypes::adolc_taped>(); + test_scalar<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_01_1.output b/tests/adolc/helper_scalar_single_component_01_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_01_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_01_2.cc b/tests/adolc/helper_scalar_single_component_01_2.cc new file mode 100644 index 0000000000..73db2465f4 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_01_2.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (scalar) system using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_single_component_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_scalar<2, double, AD::NumberTypes::adolc_tapeless>(); + test_scalar<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_01_2.output b/tests/adolc/helper_scalar_single_component_01_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_01_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_02_1.cc b/tests/adolc/helper_scalar_single_component_02_1.cc new file mode 100644 index 0000000000..c25c0b512f --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_02_1.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (vector) system using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_single_component_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_vector<2, double, AD::NumberTypes::adolc_taped>(); + test_vector<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_02_1.output b/tests/adolc/helper_scalar_single_component_02_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_02_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_02_2.cc b/tests/adolc/helper_scalar_single_component_02_2.cc new file mode 100644 index 0000000000..ef10d022ed --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_02_2.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (vector) system using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_single_component_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_vector<2, double, AD::NumberTypes::adolc_tapeless>(); + test_vector<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_02_2.output b/tests/adolc/helper_scalar_single_component_02_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_02_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_03_1.cc b/tests/adolc/helper_scalar_single_component_03_1.cc new file mode 100644 index 0000000000..981d45207e --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_03_1.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (tensor) system using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_single_component_03.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor<2, double, AD::NumberTypes::adolc_taped>(); + test_tensor<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_03_1.output b/tests/adolc/helper_scalar_single_component_03_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_03_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_03_2.cc b/tests/adolc/helper_scalar_single_component_03_2.cc new file mode 100644 index 0000000000..3118095809 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_03_2.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (tensor) system using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_single_component_03.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor<2, double, AD::NumberTypes::adolc_tapeless>(); + test_tensor<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_03_2.output b/tests/adolc/helper_scalar_single_component_03_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_03_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_04_1.cc b/tests/adolc/helper_scalar_single_component_04_1.cc new file mode 100644 index 0000000000..999ca53df3 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_04_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (symmetric tensor) system using a helper +// class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_scalar_single_component_04.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_04_1.output b/tests/adolc/helper_scalar_single_component_04_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_04_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_scalar_single_component_04_2.cc b/tests/adolc/helper_scalar_single_component_04_2.cc new file mode 100644 index 0000000000..feebc1cacf --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_04_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (symmetric tensor) system using a helper +// class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_scalar_single_component_04.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_scalar_single_component_04_2.output b/tests/adolc/helper_scalar_single_component_04_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_scalar_single_component_04_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_tape_index_01.cc b/tests/adolc/helper_tape_index_01.cc new file mode 100644 index 0000000000..6245a904cd --- /dev/null +++ b/tests/adolc/helper_tape_index_01.cc @@ -0,0 +1,338 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Verify that the maximum number of indices that can be used to record tapes is +// indeed that stored within the HelperBase class +// +// AD number type: ADOL-C taped + +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; +namespace AD = dealii::Differentiation::AD; + +// Function and its derivatives +template +struct FunctionsTestScalarScalarCoupled +{ + static NumberType + psi(const NumberType &s1, const NumberType &s2) + { + return 2.0 * std::pow(s1, 4) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds1(const NumberType &s1, const NumberType &s2) + { + return 8.0 * std::pow(s1, 3) * std::pow(s2, 3); + }; + + static NumberType + dpsi_ds2(const NumberType &s1, const NumberType &s2) + { + return 6.0 * std::pow(s1, 4) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 2) * std::pow(s2, 3); + }; + + static NumberType + d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2) + { + return 24.0 * std::pow(s1, 3) * std::pow(s2, 2); + }; + + static NumberType + d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2) + { + return d2psi_ds2_ds1(s1, s2); + }; + + static NumberType + d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2) + { + return 12.0 * std::pow(s1, 4) * std::pow(s2, 1); + }; +}; + +template +void +test_scalar_scalar_coupled_one_run(const number_t s1, + const number_t s2, + const tape_index_t tape_no) +{ + typedef AD::ScalarFunction ADHelper; + typedef typename ADHelper::ad_type ADNumberType; + typedef typename ADHelper::scalar_type ScalarNumberType; + + std::cout + // << "*** Test variables: Scalar + Scalar (coupled), " + // << "dim = " << Utilities::to_string(dim) << ", " + // << "Type code: " << static_cast(ad_type_code) << ". " + << "Tape number: " << tape_no << std::endl; + + // Values computed from the AD energy function + ScalarNumberType psi; + Vector Dpsi; + FullMatrix D2psi; + + // Function and its derivatives + typedef FunctionsTestScalarScalarCoupled func_ad; + + // Setup the variable components and choose a value at which to + // evaluate the tape + const FEValuesExtractors::Scalar s1_dof(0); + const FEValuesExtractors::Scalar s2_dof(1); + const unsigned int n_AD_components = 2; + ADHelper ad_helper(n_AD_components); + ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the + // default values + const bool write_tapes_to_file = true; // This flag will help check whether + // tape indices are universally valid + + // Configure tape + const bool is_recording = + ad_helper.start_recording_operations(tape_no, + false /*overwrite_tape*/, + true /*keep*/); + if (is_recording == true) + { + ad_helper.register_independent_variable(s1, s1_dof); + ad_helper.register_independent_variable(s2, s2_dof); + + const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof); + const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof); + + const ADNumberType psi(func_ad::psi(s1_ad, s2_ad)); + + ad_helper.register_dependent_variable(psi); + ad_helper.stop_recording_operations(write_tapes_to_file); + + // std::cout << "Taped data..." << std::endl; + // std::cout << "independent variable values: " << std::flush; + // ad_helper.print_values(std::cout); + // std::cout << "s1_ad: " << s1_ad << std::endl; + // std::cout << "s2_ad: " << s2_ad << std::endl; + // std::cout << "psi: " << psi << std::endl; + // std::cout << std::endl; + } + else + { + // std::cout << "Using tape with different values for independent + // variables..." << std::endl; + // Set a new evaluation point + ad_helper.activate_recorded_tape(tape_no); + ad_helper.set_independent_variable(s1, s1_dof); + ad_helper.set_independent_variable(s2, s2_dof); + } + + // std::cout << "independent variable values: " << std::flush; + // ad_helper.print_values(std::cout); + + // Compute the function value, gradient and hessian for the new evaluation + // point + psi = ad_helper.compute_value(); + ad_helper.compute_gradient(Dpsi); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + ad_helper.compute_hessian(D2psi); + } + + // Output the full stored function, gradient vector and hessian matrix + // std::cout << "psi: " << psi << std::endl; + // std::cout << "Dpsi: \n"; + // Dpsi.print(std::cout); + // if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + // { + // std::cout << "D2psi: \n"; + // D2psi.print_formatted(std::cout,3,true,0,"0.0"); + // } + + // Extract components of the solution + const ScalarNumberType dpsi_ds1 = + ad_helper.extract_gradient_component(Dpsi, s1_dof); + const ScalarNumberType dpsi_ds2 = + ad_helper.extract_gradient_component(Dpsi, s2_dof); + // std::cout + // << "extracted Dpsi (s1): " << dpsi_ds1 << "\n" + // << "extracted Dpsi (s2): " << dpsi_ds2 << "\n"; + + // Verify the result + typedef FunctionsTestScalarScalarCoupled func; + static const ScalarNumberType tol = + 1e8 * std::numeric_limits::epsilon(); + // std::cout << "dpsi_ds1: " << dpsi_ds1 << std::endl; + // std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1,s2) << + // std::endl; std::cout << "diff: " << std::abs(dpsi_ds1 - + // func::dpsi_ds1(s1,s2)) << std::endl; std::cout << "dpsi_ds2: " + // << dpsi_ds2 << std::endl; std::cout << "func::dpsi_ds2(s1,s2): " << + // func::dpsi_ds2(s1,s2) << std::endl; std::cout << "diff: " << + // std::abs(dpsi_ds2 - func::dpsi_ds2(s1,s2)) << std::endl; + Assert(std::abs(psi - func::psi(s1, s2)) < tol, + ExcMessage("No match for function value.")); + Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol, + ExcMessage("No match for first derivative.")); + if (AD::ADNumberTraits::n_supported_derivative_levels >= 2) + { + const ScalarNumberType d2psi_ds1_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds1 = + ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof); + const ScalarNumberType d2psi_ds1_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof); + const ScalarNumberType d2psi_ds2_ds2 = + ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof); + // std::cout + // << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n" + // << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n" + // << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n" + // << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n" + // << std::endl; + // std::cout << "d2psi_ds1_ds1: " << d2psi_ds1_ds1 << + // std::endl; std::cout << "func::d2psi_ds1_ds1(s1,s2): " << + // func::d2psi_ds1_ds1(s1,s2) << std::endl; std::cout << "diff: " << + // std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1,s2)) << std::endl; + // std::cout << "d2psi_ds2_ds1: " << d2psi_ds2_ds1 << + // std::endl; std::cout << "func::d2psi_ds2_ds1(s1,s2): " << + // func::d2psi_ds2_ds1(s1,s2) << std::endl; std::cout << "diff: " << + // std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1,s2)) << std::endl; + // std::cout << "d2psi_ds1_ds2: " << d2psi_ds1_ds2 << + // std::endl; std::cout << "func::d2psi_ds1_ds2(s1,s2): " << + // func::d2psi_ds1_ds2(s1,s2) << std::endl; std::cout << "diff: " << + // std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1,s2)) << std::endl; + // std::cout << "d2psi_ds2_ds2: " << d2psi_ds2_ds2 << + // std::endl; std::cout << "func::d2psi_ds2_ds2(s1,s2): " << + // func::d2psi_ds2_ds2(s1,s2) << std::endl; std::cout << "diff: " << + // std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1,s2)) << std::endl; + Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol, + ExcMessage("No match for second derivative.")); + } +} + + + +template +void +test_scalar_scalar_coupled(const tape_index_t tape_index_begin, + const tape_index_t tape_index_end) +{ + for (tape_index_t tape_no = tape_index_begin; tape_no < tape_index_end; + ++tape_no) + { + const number_t s1 = 3.1 * (1.0 + tape_no / 10.0); + const number_t s2 = 5.9 * (1.0 + tape_no / 10.0); + test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>(s1, + s2, + tape_no); + } +} + + +int +main() +{ + initlog(); + + using number_t = double; + const AD::NumberTypes ad_type_code = AD::NumberTypes::adolc_taped; + using ADNumberType = typename AD::HelperBase::ad_type; + using tape_index_t = typename AD::Types::tape_index; + + const tape_index_t invalid_tape_index = + AD::Numbers::invalid_tape_index; + const tape_index_t max_tape_index = AD::Numbers::max_tape_index; + + deallog.push("Valid range"); + { + test_scalar_scalar_coupled<2, number_t, ad_type_code>( + static_cast(invalid_tape_index + 1), max_tape_index); + test_scalar_scalar_coupled<3, number_t, ad_type_code>( + static_cast(invalid_tape_index + 1), max_tape_index); + deallog << "OK" << std::endl; + } + deallog.pop(); + + // We expect the following two checks to fail + deal_II_exceptions::disable_abort_on_exception(); + + deallog.push("Invalid index"); + { + try + { + test_scalar_scalar_coupled<2, number_t, ad_type_code>( + invalid_tape_index, + static_cast(invalid_tape_index + 1)); + test_scalar_scalar_coupled<3, number_t, ad_type_code>( + invalid_tape_index, + static_cast(invalid_tape_index + 1)); + deallog << "Expected failure" << std::endl; + } + catch (...) + { + deallog << "OK" << std::endl; + } + } + deallog.pop(); + + deallog.push("Max index"); + { + try + { + test_scalar_scalar_coupled<2, number_t, ad_type_code>( + max_tape_index, static_cast(max_tape_index + 1)); + test_scalar_scalar_coupled<3, number_t, ad_type_code>( + max_tape_index, static_cast(max_tape_index + 1)); + deallog << "Expected failure" << std::endl; + } + catch (...) + { + deallog << "OK" << std::endl; + } + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_tape_index_01.output b/tests/adolc/helper_tape_index_01.output new file mode 100644 index 0000000000..666b4bc9e5 --- /dev/null +++ b/tests/adolc/helper_tape_index_01.output @@ -0,0 +1,5 @@ + +DEAL:Valid range::OK +DEAL:Invalid index::OK +DEAL:Max index::OK +DEAL::OK diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc new file mode 100644 index 0000000000..687ea6b3fb --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (scalar) system using a helper class +// This is the equivalent of helper_scalar_single_component_01, but for +// the other (vector) helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_vector_1_indep_1_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc new file mode 100644 index 0000000000..525b2c7c37 --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a single component (scalar) system using a helper class +// This is the equivalent of helper_scalar_single_component_01, but for +// the other (vector) helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_vector_1_indep_1_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc new file mode 100644 index 0000000000..b63f18e1d8 --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a vector of 2 dependent and 1 independent variables +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_vector_1_indep_2_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc new file mode 100644 index 0000000000..da1461400d --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a vector of 2 dependent and 1 independent variables +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_vector_1_indep_2_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc new file mode 100644 index 0000000000..a725273f70 --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a vector of 1 dependent and 2 independent variables +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_vector_2_indep_1_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc new file mode 100644 index 0000000000..33713807f4 --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a vector of 1 dependent and 2 independent variables +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_vector_2_indep_1_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc new file mode 100644 index 0000000000..264f6c28fa --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a vector of 2 dependent and 2 independent variables +// using a helper class +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/helper_vector_2_indep_2_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc new file mode 100644 index 0000000000..969a2faa5b --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Evaluation of a vector of 2 dependent and 2 independent variables +// using a helper class +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/helper_vector_2_indep_2_dep_vars_01.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>(); + test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/step-44-helper_scalar_01_1.cc b/tests/adolc/step-44-helper_scalar_01_1.cc new file mode 100644 index 0000000000..b1de2d10b3 --- /dev/null +++ b/tests/adolc/step-44-helper_scalar_01_1.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// This is a modified version of step-44, which tests the implementation of +// QP-level auto-differentiation. +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/step-44-helper_scalar_01.h" +#include "../tests.h" + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + namespace AD = dealii::Differentiation::AD; + using namespace Step44; + try + { + const unsigned int dim = 2; + Solid solid( + SOURCE_DIR "/prm/parameters-step-44.prm"); + solid.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +} diff --git a/tests/adolc/step-44-helper_scalar_01_1.output b/tests/adolc/step-44-helper_scalar_01_1.output new file mode 100644 index 0000000000..8ba79aecc8 --- /dev/null +++ b/tests/adolc/step-44-helper_scalar_01_1.output @@ -0,0 +1,11 @@ + +DEAL::Timestep 1: 0.00 -0.000119 +DEAL::Timestep 2: 0.00 -0.000228 +DEAL::Timestep 3: 0.00 -0.000322 +DEAL::Timestep 4: 0.00 -0.000402 +DEAL::Timestep 5: 0.00 -0.000470 +DEAL::Timestep 6: 0.00 -0.000527 +DEAL::Timestep 7: 0.00 -0.000576 +DEAL::Timestep 8: 0.00 -0.000617 +DEAL::Timestep 9: 0.00 -0.000653 +DEAL::Timestep 10: 0.00 -0.000683 diff --git a/tests/adolc/step-44-helper_vector_01_1.cc b/tests/adolc/step-44-helper_vector_01_1.cc new file mode 100644 index 0000000000..27588383dd --- /dev/null +++ b/tests/adolc/step-44-helper_vector_01_1.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// This is a modified version of step-44, which tests the implementation of +// cell-level auto-differentiation (linearisation of a residual vector). +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/step-44-helper_vector_01.h" +#include "../tests.h" + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + namespace AD = dealii::Differentiation::AD; + using namespace Step44; + try + { + const unsigned int dim = 2; + Solid solid( + SOURCE_DIR "/prm/parameters-step-44.prm"); + solid.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +} diff --git a/tests/adolc/step-44-helper_vector_01_1.output b/tests/adolc/step-44-helper_vector_01_1.output new file mode 100644 index 0000000000..8ba79aecc8 --- /dev/null +++ b/tests/adolc/step-44-helper_vector_01_1.output @@ -0,0 +1,11 @@ + +DEAL::Timestep 1: 0.00 -0.000119 +DEAL::Timestep 2: 0.00 -0.000228 +DEAL::Timestep 3: 0.00 -0.000322 +DEAL::Timestep 4: 0.00 -0.000402 +DEAL::Timestep 5: 0.00 -0.000470 +DEAL::Timestep 6: 0.00 -0.000527 +DEAL::Timestep 7: 0.00 -0.000576 +DEAL::Timestep 8: 0.00 -0.000617 +DEAL::Timestep 9: 0.00 -0.000653 +DEAL::Timestep 10: 0.00 -0.000683 diff --git a/tests/adolc/step-44-helper_vector_01_2.cc b/tests/adolc/step-44-helper_vector_01_2.cc new file mode 100644 index 0000000000..55f51f4f03 --- /dev/null +++ b/tests/adolc/step-44-helper_vector_01_2.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// This is a modified version of step-44, which tests the implementation of +// cell-level auto-differentiation (linearisation of a residual vector). +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/step-44-helper_vector_01.h" +#include "../tests.h" + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + namespace AD = dealii::Differentiation::AD; + using namespace Step44; + try + { + const unsigned int dim = 2; + Solid solid( + SOURCE_DIR "/prm/parameters-step-44.prm"); + solid.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +} diff --git a/tests/adolc/step-44-helper_vector_01_2.output b/tests/adolc/step-44-helper_vector_01_2.output new file mode 100644 index 0000000000..8ba79aecc8 --- /dev/null +++ b/tests/adolc/step-44-helper_vector_01_2.output @@ -0,0 +1,11 @@ + +DEAL::Timestep 1: 0.00 -0.000119 +DEAL::Timestep 2: 0.00 -0.000228 +DEAL::Timestep 3: 0.00 -0.000322 +DEAL::Timestep 4: 0.00 -0.000402 +DEAL::Timestep 5: 0.00 -0.000470 +DEAL::Timestep 6: 0.00 -0.000527 +DEAL::Timestep 7: 0.00 -0.000576 +DEAL::Timestep 8: 0.00 -0.000617 +DEAL::Timestep 9: 0.00 -0.000653 +DEAL::Timestep 10: 0.00 -0.000683 diff --git a/tests/adolc/step-44.cc b/tests/adolc/step-44.cc new file mode 100644 index 0000000000..7e7c78f18d --- /dev/null +++ b/tests/adolc/step-44.cc @@ -0,0 +1,2128 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// This is a copy of step-44 (git rev 3f7e617) to use as a base-line for +// results produced via different approaches to be compared to. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" +namespace Step44 +{ + using namespace dealii; + namespace Parameters + { + struct FESystem + { + unsigned int poly_degree; + unsigned int quad_order; + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + void + FESystem::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + prm.declare_entry("Polynomial degree", + "2", + Patterns::Integer(0), + "Displacement system polynomial order"); + prm.declare_entry("Quadrature order", + "3", + Patterns::Integer(0), + "Gauss quadrature order"); + } + prm.leave_subsection(); + } + void + FESystem::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + poly_degree = prm.get_integer("Polynomial degree"); + quad_order = prm.get_integer("Quadrature order"); + } + prm.leave_subsection(); + } + struct Geometry + { + unsigned int global_refinement; + double scale; + double p_p0; + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + void + Geometry::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + prm.declare_entry("Global refinement", + "2", + Patterns::Integer(0), + "Global refinement level"); + prm.declare_entry("Grid scale", + "1e-3", + Patterns::Double(0.0), + "Global grid scaling factor"); + prm.declare_entry("Pressure ratio p/p0", + "100", + Patterns::Selection("20|40|60|80|100"), + "Ratio of applied pressure to reference pressure"); + } + prm.leave_subsection(); + } + void + Geometry::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + global_refinement = prm.get_integer("Global refinement"); + scale = prm.get_double("Grid scale"); + p_p0 = prm.get_double("Pressure ratio p/p0"); + } + prm.leave_subsection(); + } + struct Materials + { + double nu; + double mu; + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + void + Materials::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + prm.declare_entry("Poisson's ratio", + "0.4999", + Patterns::Double(-1.0, 0.5), + "Poisson's ratio"); + prm.declare_entry("Shear modulus", + "80.194e6", + Patterns::Double(), + "Shear modulus"); + } + prm.leave_subsection(); + } + void + Materials::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + nu = prm.get_double("Poisson's ratio"); + mu = prm.get_double("Shear modulus"); + } + prm.leave_subsection(); + } + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + bool use_static_condensation; + std::string preconditioner_type; + double preconditioner_relaxation; + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + void + LinearSolver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + prm.declare_entry("Solver type", + "CG", + Patterns::Selection("CG|Direct"), + "Type of solver used to solve the linear system"); + prm.declare_entry("Residual", + "1e-6", + Patterns::Double(0.0), + "Linear solver residual (scaled by residual norm)"); + prm.declare_entry( + "Max iteration multiplier", + "1", + Patterns::Double(0.0), + "Linear solver iterations (multiples of the system matrix size)"); + prm.declare_entry("Use static condensation", + "true", + Patterns::Bool(), + "Solve the full block system or a reduced problem"); + prm.declare_entry("Preconditioner type", + "ssor", + Patterns::Selection("jacobi|ssor"), + "Type of preconditioner"); + prm.declare_entry("Preconditioner relaxation", + "0.65", + Patterns::Double(0.0), + "Preconditioner relaxation value"); + } + prm.leave_subsection(); + } + void + LinearSolver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + type_lin = prm.get("Solver type"); + tol_lin = prm.get_double("Residual"); + max_iterations_lin = prm.get_double("Max iteration multiplier"); + use_static_condensation = prm.get_bool("Use static condensation"); + preconditioner_type = prm.get("Preconditioner type"); + preconditioner_relaxation = prm.get_double("Preconditioner relaxation"); + } + prm.leave_subsection(); + } + struct NonlinearSolver + { + unsigned int max_iterations_NR; + double tol_f; + double tol_u; + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + void + NonlinearSolver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Nonlinear solver"); + { + prm.declare_entry("Max iterations Newton-Raphson", + "10", + Patterns::Integer(0), + "Number of Newton-Raphson iterations allowed"); + prm.declare_entry("Tolerance force", + "1.0e-9", + Patterns::Double(0.0), + "Force residual tolerance"); + prm.declare_entry("Tolerance displacement", + "1.0e-6", + Patterns::Double(0.0), + "Displacement error tolerance"); + } + prm.leave_subsection(); + } + void + NonlinearSolver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Nonlinear solver"); + { + max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); + tol_f = prm.get_double("Tolerance force"); + tol_u = prm.get_double("Tolerance displacement"); + } + prm.leave_subsection(); + } + struct Time + { + double delta_t; + double end_time; + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + void + Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "1", Patterns::Double(), "End time"); + prm.declare_entry("Time step size", + "0.1", + Patterns::Double(), + "Time step size"); + } + prm.leave_subsection(); + } + void + Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); + } + prm.leave_subsection(); + } + struct AllParameters : public FESystem, + public Geometry, + public Materials, + public LinearSolver, + public NonlinearSolver, + public Time + { + AllParameters(const std::string &input_file); + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.parse_input(input_file); + parse_parameters(prm); + } + void + AllParameters::declare_parameters(ParameterHandler &prm) + { + FESystem::declare_parameters(prm); + Geometry::declare_parameters(prm); + Materials::declare_parameters(prm); + LinearSolver::declare_parameters(prm); + NonlinearSolver::declare_parameters(prm); + Time::declare_parameters(prm); + } + void + AllParameters::parse_parameters(ParameterHandler &prm) + { + FESystem::parse_parameters(prm); + Geometry::parse_parameters(prm); + Materials::parse_parameters(prm); + LinearSolver::parse_parameters(prm); + NonlinearSolver::parse_parameters(prm); + Time::parse_parameters(prm); + } + } // namespace Parameters + template + class StandardTensors + { + public: + static const SymmetricTensor<2, dim> I; + static const SymmetricTensor<4, dim> IxI; + static const SymmetricTensor<4, dim> II; + static const SymmetricTensor<4, dim> dev_P; + }; + template + const SymmetricTensor<2, dim> + StandardTensors::I = unit_symmetric_tensor(); + template + const SymmetricTensor<4, dim> StandardTensors::IxI = outer_product(I, I); + template + const SymmetricTensor<4, dim> + StandardTensors::II = identity_tensor(); + template + const SymmetricTensor<4, dim> + StandardTensors::dev_P = deviator_tensor(); + class Time + { + public: + Time(const double time_end, const double delta_t) + : timestep(0) + , time_current(0.0) + , time_end(time_end) + , delta_t(delta_t) + {} + virtual ~Time() + {} + double + current() const + { + return time_current; + } + double + end() const + { + return time_end; + } + double + get_delta_t() const + { + return delta_t; + } + unsigned int + get_timestep() const + { + return timestep; + } + void + increment() + { + time_current += delta_t; + ++timestep; + } + + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; + }; + template + class Material_Compressible_Neo_Hook_Three_Field + { + public: + Material_Compressible_Neo_Hook_Three_Field(const double mu, const double nu) + : kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))) + , c_1(mu / 2.0) + , det_F(1.0) + , p_tilde(0.0) + , J_tilde(1.0) + , b_bar(StandardTensors::I) + { + Assert(kappa > 0, ExcInternalError()); + } + ~Material_Compressible_Neo_Hook_Three_Field() + {} + void + update_material_data(const Tensor<2, dim> &F, + const double p_tilde_in, + const double J_tilde_in) + { + det_F = determinant(F); + b_bar = std::pow(det_F, -2.0 / dim) * symmetrize(F * transpose(F)); + p_tilde = p_tilde_in; + J_tilde = J_tilde_in; + Assert(det_F > 0, ExcInternalError()); + } + SymmetricTensor<2, dim> + get_tau() + { + return get_tau_iso() + get_tau_vol(); + } + SymmetricTensor<4, dim> + get_Jc() const + { + return get_Jc_vol() + get_Jc_iso(); + } + double + get_dPsi_vol_dJ() const + { + return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde); + } + double + get_d2Psi_vol_dJ2() const + { + return ((kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde))); + } + double + get_det_F() const + { + return det_F; + } + double + get_p_tilde() const + { + return p_tilde; + } + double + get_J_tilde() const + { + return J_tilde; + } + + protected: + const double kappa; + const double c_1; + double det_F; + double p_tilde; + double J_tilde; + SymmetricTensor<2, dim> b_bar; + SymmetricTensor<2, dim> + get_tau_vol() const + { + return p_tilde * det_F * StandardTensors::I; + } + SymmetricTensor<2, dim> + get_tau_iso() const + { + return StandardTensors::dev_P * get_tau_bar(); + } + SymmetricTensor<2, dim> + get_tau_bar() const + { + return 2.0 * c_1 * b_bar; + } + SymmetricTensor<4, dim> + get_Jc_vol() const + { + return p_tilde * det_F * + (StandardTensors::IxI - (2.0 * StandardTensors::II)); + } + SymmetricTensor<4, dim> + get_Jc_iso() const + { + const SymmetricTensor<2, dim> tau_bar = get_tau_bar(); + const SymmetricTensor<2, dim> tau_iso = get_tau_iso(); + const SymmetricTensor<4, dim> tau_iso_x_I = + outer_product(tau_iso, StandardTensors::I); + const SymmetricTensor<4, dim> I_x_tau_iso = + outer_product(StandardTensors::I, tau_iso); + const SymmetricTensor<4, dim> c_bar = get_c_bar(); + return (2.0 / dim) * trace(tau_bar) * StandardTensors::dev_P - + (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso) + + StandardTensors::dev_P * c_bar * StandardTensors::dev_P; + } + SymmetricTensor<4, dim> + get_c_bar() const + { + return SymmetricTensor<4, dim>(); + } + }; + template + class PointHistory + { + public: + PointHistory() + : F_inv(StandardTensors::I) + , tau(SymmetricTensor<2, dim>()) + , d2Psi_vol_dJ2(0.0) + , dPsi_vol_dJ(0.0) + , Jc(SymmetricTensor<4, dim>()) + {} + virtual ~PointHistory() + {} + void + setup_lqp(const Parameters::AllParameters ¶meters) + { + material.reset( + new Material_Compressible_Neo_Hook_Three_Field(parameters.mu, + parameters.nu)); + update_values(Tensor<2, dim>(), 0.0, 1.0); + } + void + update_values(const Tensor<2, dim> &Grad_u_n, + const double p_tilde, + const double J_tilde) + { + const Tensor<2, dim> F = + (Tensor<2, dim>(StandardTensors::I) + Grad_u_n); + material->update_material_data(F, p_tilde, J_tilde); + F_inv = invert(F); + tau = material->get_tau(); + Jc = material->get_Jc(); + dPsi_vol_dJ = material->get_dPsi_vol_dJ(); + d2Psi_vol_dJ2 = material->get_d2Psi_vol_dJ2(); + } + double + get_J_tilde() const + { + return material->get_J_tilde(); + } + double + get_det_F() const + { + return material->get_det_F(); + } + const Tensor<2, dim> & + get_F_inv() const + { + return F_inv; + } + double + get_p_tilde() const + { + return material->get_p_tilde(); + } + const SymmetricTensor<2, dim> & + get_tau() const + { + return tau; + } + double + get_dPsi_vol_dJ() const + { + return dPsi_vol_dJ; + } + double + get_d2Psi_vol_dJ2() const + { + return d2Psi_vol_dJ2; + } + const SymmetricTensor<4, dim> & + get_Jc() const + { + return Jc; + } + + private: + std::shared_ptr> material; + Tensor<2, dim> F_inv; + SymmetricTensor<2, dim> tau; + double d2Psi_vol_dJ2; + double dPsi_vol_dJ; + SymmetricTensor<4, dim> Jc; + }; + template + class Solid + { + public: + Solid(const std::string &input_file); + virtual ~Solid(); + void + run(); + + private: + struct PerTaskData_K; + struct ScratchData_K; + struct PerTaskData_RHS; + struct ScratchData_RHS; + struct PerTaskData_SC; + struct ScratchData_SC; + struct PerTaskData_UQPH; + struct ScratchData_UQPH; + void + make_grid(); + void + system_setup(); + void + determine_component_extractors(); + void + assemble_system_tangent(); + void + assemble_system_tangent_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_K & scratch, + PerTaskData_K & data) const; + void + copy_local_to_global_K(const PerTaskData_K &data); + void + assemble_system_rhs(); + void + assemble_system_rhs_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_RHS & scratch, + PerTaskData_RHS & data) const; + void + copy_local_to_global_rhs(const PerTaskData_RHS &data); + void + assemble_sc(); + void + assemble_sc_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_SC & scratch, + PerTaskData_SC & data); + void + copy_local_to_global_sc(const PerTaskData_SC &data); + void + make_constraints(const int &it_nr); + void + setup_qph(); + void + update_qph_incremental(const BlockVector &solution_delta); + void + update_qph_incremental_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & data); + void + copy_local_to_global_UQPH(const PerTaskData_UQPH & /*data*/) + {} + void + solve_nonlinear_timestep(BlockVector &solution_delta); + std::pair + solve_linear_system(BlockVector &newton_update); + BlockVector + get_total_solution(const BlockVector &solution_delta) const; + void + output_results() const; + Parameters::AllParameters parameters; + double vol_reference; + Triangulation triangulation; + Time time; + mutable TimerOutput timer; + CellDataStorage::cell_iterator, + PointHistory> + quadrature_point_history; + const unsigned int degree; + const FESystem fe; + DoFHandler dof_handler_ref; + const unsigned int dofs_per_cell; + const FEValuesExtractors::Vector u_fe; + const FEValuesExtractors::Scalar p_fe; + const FEValuesExtractors::Scalar J_fe; + static const unsigned int n_blocks = 3; + static const unsigned int n_components = dim + 2; + static const unsigned int first_u_component = 0; + static const unsigned int p_component = dim; + static const unsigned int J_component = dim + 1; + enum + { + u_dof = 0, + p_dof = 1, + J_dof = 2 + }; + std::vector dofs_per_block; + std::vector element_indices_u; + std::vector element_indices_p; + std::vector element_indices_J; + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const unsigned int n_q_points_f; + ConstraintMatrix constraints; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix tangent_matrix; + BlockVector system_rhs; + BlockVector solution_n; + struct Errors + { + Errors() + : norm(1.0) + , u(1.0) + , p(1.0) + , J(1.0) + {} + void + reset() + { + norm = 1.0; + u = 1.0; + p = 1.0; + J = 1.0; + } + void + normalise(const Errors &rhs) + { + if (rhs.norm != 0.0) + norm /= rhs.norm; + if (rhs.u != 0.0) + u /= rhs.u; + if (rhs.p != 0.0) + p /= rhs.p; + if (rhs.J != 0.0) + J /= rhs.J; + } + double norm, u, p, J; + }; + Errors error_residual, error_residual_0, error_residual_norm, error_update, + error_update_0, error_update_norm; + void + get_error_residual(Errors &error_residual); + void + get_error_update(const BlockVector &newton_update, + Errors & error_update); + std::pair + get_error_dilation() const; + double + compute_vol_current() const; + static void + print_conv_header(); + void + print_conv_footer(); + }; + template + Solid::Solid(const std::string &input_file) + : parameters(input_file) + , triangulation(Triangulation::maximum_smoothing) + , time(parameters.end_time, parameters.delta_t) + , timer(std::cout, TimerOutput::never, TimerOutput::wall_times) + , degree(parameters.poly_degree) + , fe(FE_Q(parameters.poly_degree), + dim, // displacement + FE_DGPMonomial(parameters.poly_degree - 1), + 1, // pressure + FE_DGPMonomial(parameters.poly_degree - 1), + 1) + , // dilatation + dof_handler_ref(triangulation) + , dofs_per_cell(fe.dofs_per_cell) + , u_fe(first_u_component) + , p_fe(p_component) + , J_fe(J_component) + , dofs_per_block(n_blocks) + , qf_cell(parameters.quad_order) + , qf_face(parameters.quad_order) + , n_q_points(qf_cell.size()) + , n_q_points_f(qf_face.size()) + { + Assert(dim == 2 || dim == 3, + ExcMessage("This problem only works in 2 or 3 space dimensions.")); + determine_component_extractors(); + } + template + Solid::~Solid() + { + dof_handler_ref.clear(); + } + template + void + Solid::run() + { + make_grid(); + system_setup(); + { + ConstraintMatrix constraints; + constraints.close(); + const ComponentSelectFunction J_mask(J_component, n_components); + VectorTools::project(dof_handler_ref, + constraints, + QGauss(degree + 2), + J_mask, + solution_n); + } + // output_results(); + time.increment(); + BlockVector solution_delta(dofs_per_block); + while (time.current() < time.end()) + { + solution_delta = 0.0; + solve_nonlinear_timestep(solution_delta); + solution_n += solution_delta; + // output_results(); + // Output displacement at centre of traction surface + { + const Point soln_pt( + dim == 3 ? Point(0.0, 1.0, 0.0) * parameters.scale : + Point(0.0, 1.0) * parameters.scale); + typename DoFHandler::active_cell_iterator + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); + for (; cell != endc; ++cell) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; + ++v) + if (cell->vertex(v).distance(soln_pt) < 1e-6 * parameters.scale) + { + Tensor<1, dim> soln; + for (unsigned int d = 0; d < dim; ++d) + soln[d] = solution_n(cell->vertex_dof_index(v, u_dof + d)); + deallog << "Timestep " << time.get_timestep() << ": " << soln + << std::endl; + } + } + time.increment(); + } + } + template + struct Solid::PerTaskData_K + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + PerTaskData_K(const unsigned int dofs_per_cell) + : cell_matrix(dofs_per_cell, dofs_per_cell) + , local_dof_indices(dofs_per_cell) + {} + void + reset() + { + cell_matrix = 0.0; + } + }; + template + struct Solid::ScratchData_K + { + FEValues fe_values_ref; + std::vector> Nx; + std::vector>> grad_Nx; + std::vector>> symm_grad_Nx; + ScratchData_K(const FiniteElement &fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell) + : fe_values_ref(fe_cell, qf_cell, uf_cell) + , Nx(qf_cell.size(), std::vector(fe_cell.dofs_per_cell)) + , grad_Nx(qf_cell.size(), + std::vector>(fe_cell.dofs_per_cell)) + , symm_grad_Nx(qf_cell.size(), + std::vector>( + fe_cell.dofs_per_cell)) + {} + ScratchData_K(const ScratchData_K &rhs) + : fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + , Nx(rhs.Nx) + , grad_Nx(rhs.grad_Nx) + , symm_grad_Nx(rhs.symm_grad_Nx) + {} + void + reset() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert(grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + grad_Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + }; + template + struct Solid::PerTaskData_RHS + { + Vector cell_rhs; + std::vector local_dof_indices; + PerTaskData_RHS(const unsigned int dofs_per_cell) + : cell_rhs(dofs_per_cell) + , local_dof_indices(dofs_per_cell) + {} + void + reset() + { + cell_rhs = 0.0; + } + }; + template + struct Solid::ScratchData_RHS + { + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; + std::vector> Nx; + std::vector>> symm_grad_Nx; + ScratchData_RHS(const FiniteElement &fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const QGauss & qf_face, + const UpdateFlags uf_face) + : fe_values_ref(fe_cell, qf_cell, uf_cell) + , fe_face_values_ref(fe_cell, qf_face, uf_face) + , Nx(qf_cell.size(), std::vector(fe_cell.dofs_per_cell)) + , symm_grad_Nx(qf_cell.size(), + std::vector>( + fe_cell.dofs_per_cell)) + {} + ScratchData_RHS(const ScratchData_RHS &rhs) + : fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + , fe_face_values_ref(rhs.fe_face_values_ref.get_fe(), + rhs.fe_face_values_ref.get_quadrature(), + rhs.fe_face_values_ref.get_update_flags()) + , Nx(rhs.Nx) + , symm_grad_Nx(rhs.symm_grad_Nx) + {} + void + reset() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + }; + template + struct Solid::PerTaskData_SC + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + FullMatrix k_orig; + FullMatrix k_pu; + FullMatrix k_pJ; + FullMatrix k_JJ; + FullMatrix k_pJ_inv; + FullMatrix k_bbar; + FullMatrix A; + FullMatrix B; + FullMatrix C; + PerTaskData_SC(const unsigned int dofs_per_cell, + const unsigned int n_u, + const unsigned int n_p, + const unsigned int n_J) + : cell_matrix(dofs_per_cell, dofs_per_cell) + , local_dof_indices(dofs_per_cell) + , k_orig(dofs_per_cell, dofs_per_cell) + , k_pu(n_p, n_u) + , k_pJ(n_p, n_J) + , k_JJ(n_J, n_J) + , k_pJ_inv(n_p, n_J) + , k_bbar(n_u, n_u) + , A(n_J, n_u) + , B(n_J, n_u) + , C(n_p, n_u) + {} + void + reset() + {} + }; + template + struct Solid::ScratchData_SC + { + void + reset() + {} + }; + template + struct Solid::PerTaskData_UQPH + { + void + reset() + {} + }; + template + struct Solid::ScratchData_UQPH + { + const BlockVector & solution_total; + std::vector> solution_grads_u_total; + std::vector solution_values_p_total; + std::vector solution_values_J_total; + FEValues fe_values_ref; + ScratchData_UQPH(const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const BlockVector &solution_total) + : solution_total(solution_total) + , solution_grads_u_total(qf_cell.size()) + , solution_values_p_total(qf_cell.size()) + , solution_values_J_total(qf_cell.size()) + , fe_values_ref(fe_cell, qf_cell, uf_cell) + {} + ScratchData_UQPH(const ScratchData_UQPH &rhs) + : solution_total(rhs.solution_total) + , solution_grads_u_total(rhs.solution_grads_u_total) + , solution_values_p_total(rhs.solution_values_p_total) + , solution_values_J_total(rhs.solution_values_J_total) + , fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + {} + void + reset() + { + const unsigned int n_q_points = solution_grads_u_total.size(); + for (unsigned int q = 0; q < n_q_points; ++q) + { + solution_grads_u_total[q] = 0.0; + solution_values_p_total[q] = 0.0; + solution_values_J_total[q] = 0.0; + } + } + }; + template + void + Solid::make_grid() + { + GridGenerator::hyper_rectangle( + triangulation, + (dim == 3 ? Point(0.0, 0.0, 0.0) : Point(0.0, 0.0)), + (dim == 3 ? Point(1.0, 1.0, 1.0) : Point(1.0, 1.0)), + true); + GridTools::scale(parameters.scale, triangulation); + triangulation.refine_global(std::max(1U, parameters.global_refinement)); + vol_reference = GridTools::volume(triangulation); + std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl; + typename Triangulation::active_cell_iterator cell = triangulation + .begin_active(), + endc = + triangulation.end(); + for (; cell != endc; ++cell) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + if (cell->face(face)->at_boundary() == true && + cell->face(face)->center()[1] == 1.0 * parameters.scale) + { + if (dim == 3) + { + if (cell->face(face)->center()[0] < 0.5 * parameters.scale && + cell->face(face)->center()[2] < 0.5 * parameters.scale) + cell->face(face)->set_boundary_id(6); + } + else + { + if (cell->face(face)->center()[0] < 0.5 * parameters.scale) + cell->face(face)->set_boundary_id(6); + } + } + } + } + template + void + Solid::system_setup() + { + timer.enter_subsection("Setup system"); + std::vector block_component(n_components, + u_dof); // Displacement + block_component[p_component] = p_dof; // Pressure + block_component[J_component] = J_dof; // Dilatation + dof_handler_ref.distribute_dofs(fe); + DoFRenumbering::Cuthill_McKee(dof_handler_ref); + DoFRenumbering::component_wise(dof_handler_ref, block_component); + DoFTools::count_dofs_per_block(dof_handler_ref, + dofs_per_block, + block_component); + std::cout << "Triangulation:" + << "\n\t Number of active cells: " + << triangulation.n_active_cells() + << "\n\t Number of degrees of freedom: " + << dof_handler_ref.n_dofs() << std::endl; + tangent_matrix.clear(); + { + const types::global_dof_index n_dofs_u = dofs_per_block[u_dof]; + const types::global_dof_index n_dofs_p = dofs_per_block[p_dof]; + const types::global_dof_index n_dofs_J = dofs_per_block[J_dof]; + BlockDynamicSparsityPattern dsp(n_blocks, n_blocks); + dsp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u); + dsp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p); + dsp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J); + dsp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u); + dsp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p); + dsp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J); + dsp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u); + dsp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p); + dsp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J); + dsp.collect_sizes(); + Table<2, DoFTools::Coupling> coupling(n_components, n_components); + for (unsigned int ii = 0; ii < n_components; ++ii) + for (unsigned int jj = 0; jj < n_components; ++jj) + if (((ii < p_component) && (jj == J_component)) || + ((ii == J_component) && (jj < p_component)) || + ((ii == p_component) && (jj == p_component))) + coupling[ii][jj] = DoFTools::none; + else + coupling[ii][jj] = DoFTools::always; + DoFTools::make_sparsity_pattern( + dof_handler_ref, coupling, dsp, constraints, false); + sparsity_pattern.copy_from(dsp); + } + tangent_matrix.reinit(sparsity_pattern); + system_rhs.reinit(dofs_per_block); + system_rhs.collect_sizes(); + solution_n.reinit(dofs_per_block); + solution_n.collect_sizes(); + setup_qph(); + timer.leave_subsection(); + } + template + void + Solid::determine_component_extractors() + { + element_indices_u.clear(); + element_indices_p.clear(); + element_indices_J.clear(); + for (unsigned int k = 0; k < fe.dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; + if (k_group == u_dof) + element_indices_u.push_back(k); + else if (k_group == p_dof) + element_indices_p.push_back(k); + else if (k_group == J_dof) + element_indices_J.push_back(k); + else + { + Assert(k_group <= J_dof, ExcInternalError()); + } + } + } + template + void + Solid::setup_qph() + { + std::cout << " Setting up quadrature point data..." << std::endl; + quadrature_point_history.initialize(triangulation.begin_active(), + triangulation.end(), + n_q_points); + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + cell != triangulation.end(); + ++cell) + { + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->setup_lqp(parameters); + } + } + template + void + Solid::update_qph_incremental(const BlockVector &solution_delta) + { + timer.enter_subsection("Update QPH data"); + std::cout << " UQPH " << std::flush; + const BlockVector solution_total( + get_total_solution(solution_delta)); + const UpdateFlags uf_UQPH(update_values | update_gradients); + PerTaskData_UQPH per_task_data_UQPH; + ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total); + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::update_qph_incremental_one_cell, + &Solid::copy_local_to_global_UQPH, + scratch_data_UQPH, + per_task_data_UQPH); + timer.leave_subsection(); + } + template + void + Solid::update_qph_incremental_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & /*data*/) + { + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + Assert(scratch.solution_grads_u_total.size() == n_q_points, + ExcInternalError()); + Assert(scratch.solution_values_p_total.size() == n_q_points, + ExcInternalError()); + Assert(scratch.solution_values_J_total.size() == n_q_points, + ExcInternalError()); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + scratch.fe_values_ref[u_fe].get_function_gradients( + scratch.solution_total, scratch.solution_grads_u_total); + scratch.fe_values_ref[p_fe].get_function_values( + scratch.solution_total, scratch.solution_values_p_total); + scratch.fe_values_ref[J_fe].get_function_values( + scratch.solution_total, scratch.solution_values_J_total); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point], + scratch.solution_values_p_total[q_point], + scratch.solution_values_J_total[q_point]); + } + template + void + Solid::solve_nonlinear_timestep(BlockVector &solution_delta) + { + std::cout << std::endl + << "Timestep " << time.get_timestep() << " @ " << time.current() + << "s" << std::endl; + BlockVector newton_update(dofs_per_block); + error_residual.reset(); + error_residual_0.reset(); + error_residual_norm.reset(); + error_update.reset(); + error_update_0.reset(); + error_update_norm.reset(); + print_conv_header(); + unsigned int newton_iteration = 0; + for (; newton_iteration < parameters.max_iterations_NR; ++newton_iteration) + { + std::cout << " " << std::setw(2) << newton_iteration << " " + << std::flush; + tangent_matrix = 0.0; + system_rhs = 0.0; + assemble_system_rhs(); + get_error_residual(error_residual); + if (newton_iteration == 0) + error_residual_0 = error_residual; + error_residual_norm = error_residual; + error_residual_norm.normalise(error_residual_0); + if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u && + error_residual_norm.u <= parameters.tol_f) + { + std::cout << " CONVERGED! " << std::endl; + print_conv_footer(); + break; + } + assemble_system_tangent(); + make_constraints(newton_iteration); + constraints.condense(tangent_matrix, system_rhs); + const std::pair lin_solver_output = + solve_linear_system(newton_update); + get_error_update(newton_update, error_update); + if (newton_iteration == 0) + error_update_0 = error_update; + error_update_norm = error_update; + error_update_norm.normalise(error_update_0); + solution_delta += newton_update; + update_qph_incremental(solution_delta); + std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7) + << std::scientific << lin_solver_output.first << " " + << lin_solver_output.second << " " + << error_residual_norm.norm << " " << error_residual_norm.u + << " " << error_residual_norm.p << " " + << error_residual_norm.J << " " << error_update_norm.norm + << " " << error_update_norm.u << " " << error_update_norm.p + << " " << error_update_norm.J << " " << std::endl; + } + AssertThrow(newton_iteration <= parameters.max_iterations_NR, + ExcMessage("No convergence in nonlinear solver!")); + } + template + void + Solid::print_conv_header() + { + static const unsigned int l_width = 155; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + std::cout << " SOLVER STEP " + << " | LIN_IT LIN_RES RES_NORM " + << " RES_U RES_P RES_J NU_NORM " + << " NU_U NU_P NU_J " << std::endl; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + } + template + void + Solid::print_conv_footer() + { + static const unsigned int l_width = 155; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + const std::pair error_dil = get_error_dilation(); + std::cout << "Relative errors:" << std::endl + << "Displacement:\t" << error_update.u / error_update_0.u + << std::endl + << "Force: \t\t" << error_residual.u / error_residual_0.u + << std::endl + << "Dilatation:\t" << error_dil.first << std::endl + << "v / V_0:\t" << error_dil.second * vol_reference << " / " + << vol_reference << " = " << error_dil.second << std::endl; + } + template + double + Solid::compute_vol_current() const + { + double vol_current = 0.0; + FEValues fe_values_ref(fe, qf_cell, update_JxW_values); + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + cell != triangulation.end(); + ++cell) + { + fe_values_ref.reinit(cell); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double det_F_qp = lqph[q_point]->get_det_F(); + const double JxW = fe_values_ref.JxW(q_point); + vol_current += det_F_qp * JxW; + } + } + Assert(vol_current > 0.0, ExcInternalError()); + return vol_current; + } + template + std::pair + Solid::get_error_dilation() const + { + double dil_L2_error = 0.0; + FEValues fe_values_ref(fe, qf_cell, update_JxW_values); + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + cell != triangulation.end(); + ++cell) + { + fe_values_ref.reinit(cell); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double det_F_qp = lqph[q_point]->get_det_F(); + const double J_tilde_qp = lqph[q_point]->get_J_tilde(); + const double the_error_qp_squared = + std::pow((det_F_qp - J_tilde_qp), 2); + const double JxW = fe_values_ref.JxW(q_point); + dil_L2_error += the_error_qp_squared * JxW; + } + } + return std::make_pair(std::sqrt(dil_L2_error), + compute_vol_current() / vol_reference); + } + template + void + Solid::get_error_residual(Errors &error_residual) + { + BlockVector error_res(dofs_per_block); + for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i) + if (!constraints.is_constrained(i)) + error_res(i) = system_rhs(i); + error_residual.norm = error_res.l2_norm(); + error_residual.u = error_res.block(u_dof).l2_norm(); + error_residual.p = error_res.block(p_dof).l2_norm(); + error_residual.J = error_res.block(J_dof).l2_norm(); + } + template + void + Solid::get_error_update(const BlockVector &newton_update, + Errors & error_update) + { + BlockVector error_ud(dofs_per_block); + for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i) + if (!constraints.is_constrained(i)) + error_ud(i) = newton_update(i); + error_update.norm = error_ud.l2_norm(); + error_update.u = error_ud.block(u_dof).l2_norm(); + error_update.p = error_ud.block(p_dof).l2_norm(); + error_update.J = error_ud.block(J_dof).l2_norm(); + } + template + BlockVector + Solid::get_total_solution( + const BlockVector &solution_delta) const + { + BlockVector solution_total(solution_n); + solution_total += solution_delta; + return solution_total; + } + template + void + Solid::assemble_system_tangent() + { + timer.enter_subsection("Assemble tangent matrix"); + std::cout << " ASM_K " << std::flush; + tangent_matrix = 0.0; + const UpdateFlags uf_cell(update_values | update_gradients | + update_JxW_values); + PerTaskData_K per_task_data(dofs_per_cell); + ScratchData_K scratch_data(fe, qf_cell, uf_cell); + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + std::bind(&Solid::assemble_system_tangent_one_cell, + this, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3), + std::bind(&Solid::copy_local_to_global_K, + this, + std::placeholders::_1), + scratch_data, + per_task_data); + timer.leave_subsection(); + } + template + void + Solid::copy_local_to_global_K(const PerTaskData_K &data) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + tangent_matrix.add(data.local_dof_indices[i], + data.local_dof_indices[j], + data.cell_matrix(i, j)); + } + template + void + Solid::assemble_system_tangent_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_K & scratch, + PerTaskData_K & data) const + { + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv(); + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; + if (k_group == u_dof) + { + scratch.grad_Nx[q_point][k] = + scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv; + scratch.symm_grad_Nx[q_point][k] = + symmetrize(scratch.grad_Nx[q_point][k]); + } + else if (k_group == p_dof) + scratch.Nx[q_point][k] = + scratch.fe_values_ref[p_fe].value(k, q_point); + else if (k_group == J_dof) + scratch.Nx[q_point][k] = + scratch.fe_values_ref[J_fe].value(k, q_point); + else + Assert(k_group <= J_dof, ExcInternalError()); + } + } + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> tau = lqph[q_point]->get_tau(); + const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc(); + const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2(); + const double det_F = lqph[q_point]->get_det_F(); + const std::vector & N = scratch.Nx[q_point]; + const std::vector> &symm_grad_Nx = + scratch.symm_grad_Nx[q_point]; + const std::vector> &grad_Nx = scratch.grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + const unsigned int i_group = fe.system_to_base_index(i).first.first; + for (unsigned int j = 0; j <= i; ++j) + { + const unsigned int component_j = + fe.system_to_component_index(j).first; + const unsigned int j_group = + fe.system_to_base_index(j).first.first; + if ((i_group == j_group) && (i_group == u_dof)) + { + data.cell_matrix(i, j) += symm_grad_Nx[i] * + Jc // The material contribution: + * symm_grad_Nx[j] * JxW; + if (component_i == + component_j) // geometrical stress contribution + data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau * + grad_Nx[j][component_j] * JxW; + } + else if ((i_group == p_dof) && (j_group == u_dof)) + { + data.cell_matrix(i, j) += + N[i] * det_F * + (symm_grad_Nx[j] * StandardTensors::I) * JxW; + } + else if ((i_group == J_dof) && (j_group == p_dof)) + data.cell_matrix(i, j) -= N[i] * N[j] * JxW; + else if ((i_group == j_group) && (i_group == J_dof)) + data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW; + else + Assert((i_group <= J_dof) && (j_group <= J_dof), + ExcInternalError()); + } + } + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + data.cell_matrix(i, j) = data.cell_matrix(j, i); + } + template + void + Solid::assemble_system_rhs() + { + timer.enter_subsection("Assemble system right-hand side"); + std::cout << " ASM_R " << std::flush; + system_rhs = 0.0; + const UpdateFlags uf_cell(update_values | update_gradients | + update_JxW_values); + const UpdateFlags uf_face(update_values | update_normal_vectors | + update_JxW_values); + PerTaskData_RHS per_task_data(dofs_per_cell); + ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face); + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + std::bind(&Solid::assemble_system_rhs_one_cell, + this, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3), + std::bind(&Solid::copy_local_to_global_rhs, + this, + std::placeholders::_1), + scratch_data, + per_task_data); + timer.leave_subsection(); + } + template + void + Solid::copy_local_to_global_rhs(const PerTaskData_RHS &data) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i); + } + template + void + Solid::assemble_system_rhs_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_RHS & scratch, + PerTaskData_RHS & data) const + { + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv(); + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; + if (k_group == u_dof) + scratch.symm_grad_Nx[q_point][k] = symmetrize( + scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv); + else if (k_group == p_dof) + scratch.Nx[q_point][k] = + scratch.fe_values_ref[p_fe].value(k, q_point); + else if (k_group == J_dof) + scratch.Nx[q_point][k] = + scratch.fe_values_ref[J_fe].value(k, q_point); + else + Assert(k_group <= J_dof, ExcInternalError()); + } + } + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau(); + const double det_F = lqph[q_point]->get_det_F(); + const double J_tilde = lqph[q_point]->get_J_tilde(); + const double p_tilde = lqph[q_point]->get_p_tilde(); + const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ(); + const std::vector & N = scratch.Nx[q_point]; + const std::vector> &symm_grad_Nx = + scratch.symm_grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int i_group = fe.system_to_base_index(i).first.first; + if (i_group == u_dof) + data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW; + else if (i_group == p_dof) + data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW; + else if (i_group == J_dof) + data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p_tilde) * JxW; + else + Assert(i_group <= J_dof, ExcInternalError()); + } + } + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() == true && + cell->face(face)->boundary_id() == 6) + { + scratch.fe_face_values_ref.reinit(cell, face); + for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; + ++f_q_point) + { + const Tensor<1, dim> &N = + scratch.fe_face_values_ref.normal_vector(f_q_point); + static const double p0 = + -4.0 / (parameters.scale * parameters.scale); + const double time_ramp = (time.current() / time.end()); + const double pressure = p0 * parameters.p_p0 * time_ramp; + const Tensor<1, dim> traction = pressure * N; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int i_group = + fe.system_to_base_index(i).first.first; + if (i_group == u_dof) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + const double Ni = + scratch.fe_face_values_ref.shape_value(i, f_q_point); + const double JxW = + scratch.fe_face_values_ref.JxW(f_q_point); + data.cell_rhs(i) += (Ni * traction[component_i]) * JxW; + } + } + } + } + } + template + void + Solid::make_constraints(const int &it_nr) + { + std::cout << " CST " << std::flush; + if (it_nr > 1) + return; + constraints.clear(); + const bool apply_dirichlet_bc = (it_nr == 0); + const FEValuesExtractors::Scalar x_displacement(0); + const FEValuesExtractors::Scalar y_displacement(1); + { + const int boundary_id = 0; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + } + { + const int boundary_id = 2; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + } + if (dim == 3) + { + const FEValuesExtractors::Scalar z_displacement(2); + { + const int boundary_id = 3; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + } + { + const int boundary_id = 4; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + } + { + const int boundary_id = 6; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + } + } + else + { + { + const int boundary_id = 3; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + { + const int boundary_id = 6; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + } + constraints.close(); + } + template + void + Solid::assemble_sc() + { + timer.enter_subsection("Perform static condensation"); + std::cout << " ASM_SC " << std::flush; + PerTaskData_SC per_task_data(dofs_per_cell, + element_indices_u.size(), + element_indices_p.size(), + element_indices_J.size()); + ScratchData_SC scratch_data; + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::assemble_sc_one_cell, + &Solid::copy_local_to_global_sc, + scratch_data, + per_task_data); + timer.leave_subsection(); + } + template + void + Solid::copy_local_to_global_sc(const PerTaskData_SC &data) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + tangent_matrix.add(data.local_dof_indices[i], + data.local_dof_indices[j], + data.cell_matrix(i, j)); + } + template + void + Solid::assemble_sc_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_SC & scratch, + PerTaskData_SC & data) + { + data.reset(); + scratch.reset(); + cell->get_dof_indices(data.local_dof_indices); + data.k_orig.extract_submatrix_from(tangent_matrix, + data.local_dof_indices, + data.local_dof_indices); + data.k_pu.extract_submatrix_from(data.k_orig, + element_indices_p, + element_indices_u); + data.k_pJ.extract_submatrix_from(data.k_orig, + element_indices_p, + element_indices_J); + data.k_JJ.extract_submatrix_from(data.k_orig, + element_indices_J, + element_indices_J); + data.k_pJ_inv.invert(data.k_pJ); + data.k_pJ_inv.mmult(data.A, data.k_pu); + data.k_JJ.mmult(data.B, data.A); + data.k_pJ_inv.Tmmult(data.C, data.B); + data.k_pu.Tmmult(data.k_bbar, data.C); + data.k_bbar.scatter_matrix_to(element_indices_u, + element_indices_u, + data.cell_matrix); + data.k_pJ_inv.add(-1.0, data.k_pJ); + data.k_pJ_inv.scatter_matrix_to(element_indices_p, + element_indices_J, + data.cell_matrix); + } + template + std::pair + Solid::solve_linear_system(BlockVector &newton_update) + { + unsigned int lin_it = 0; + double lin_res = 0.0; + if (parameters.use_static_condensation == true) + { + BlockVector A(dofs_per_block); + BlockVector B(dofs_per_block); + { + assemble_sc(); + tangent_matrix.block(p_dof, J_dof) + .vmult(A.block(J_dof), system_rhs.block(p_dof)); + tangent_matrix.block(J_dof, J_dof) + .vmult(B.block(J_dof), A.block(J_dof)); + A.block(J_dof) = system_rhs.block(J_dof); + A.block(J_dof) -= B.block(J_dof); + tangent_matrix.block(p_dof, J_dof) + .Tvmult(A.block(p_dof), A.block(J_dof)); + tangent_matrix.block(u_dof, p_dof) + .vmult(A.block(u_dof), A.block(p_dof)); + system_rhs.block(u_dof) -= A.block(u_dof); + timer.enter_subsection("Linear solver"); + std::cout << " SLV " << std::flush; + if (parameters.type_lin == "CG") + { + const int solver_its = tangent_matrix.block(u_dof, u_dof).m() * + parameters.max_iterations_lin; + const double tol_sol = + parameters.tol_lin * system_rhs.block(u_dof).l2_norm(); + SolverControl solver_control(solver_its, tol_sol, false, false); + GrowingVectorMemory> GVM; + SolverCG> solver_CG(solver_control, GVM); + PreconditionSelector, Vector> + preconditioner(parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof)); + solver_CG.solve(tangent_matrix.block(u_dof, u_dof), + newton_update.block(u_dof), + system_rhs.block(u_dof), + preconditioner); + lin_it = solver_control.last_step(); + lin_res = solver_control.last_value(); + } + else if (parameters.type_lin == "Direct") + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(tangent_matrix.block(u_dof, u_dof)); + A_direct.vmult(newton_update.block(u_dof), + system_rhs.block(u_dof)); + lin_it = 1; + lin_res = 0.0; + } + else + Assert(false, ExcMessage("Linear solver type not implemented")); + timer.leave_subsection(); + } + constraints.distribute(newton_update); + timer.enter_subsection("Linear solver postprocessing"); + std::cout << " PP " << std::flush; + { + tangent_matrix.block(p_dof, u_dof) + .vmult(A.block(p_dof), newton_update.block(u_dof)); + A.block(p_dof) *= -1.0; + A.block(p_dof) += system_rhs.block(p_dof); + tangent_matrix.block(p_dof, J_dof) + .vmult(newton_update.block(J_dof), A.block(p_dof)); + } + constraints.distribute(newton_update); + { + tangent_matrix.block(J_dof, J_dof) + .vmult(A.block(J_dof), newton_update.block(J_dof)); + A.block(J_dof) *= -1.0; + A.block(J_dof) += system_rhs.block(J_dof); + tangent_matrix.block(p_dof, J_dof) + .Tvmult(newton_update.block(p_dof), A.block(J_dof)); + } + constraints.distribute(newton_update); + timer.leave_subsection(); + } + else + { + std::cout << " ------ " << std::flush; + timer.enter_subsection("Linear solver"); + std::cout << " SLV " << std::flush; + if (parameters.type_lin == "CG") + { + const Vector &f_u = system_rhs.block(u_dof); + const Vector &f_p = system_rhs.block(p_dof); + const Vector &f_J = system_rhs.block(J_dof); + Vector & d_u = newton_update.block(u_dof); + Vector & d_p = newton_update.block(p_dof); + Vector & d_J = newton_update.block(J_dof); + const auto K_uu = + linear_operator(tangent_matrix.block(u_dof, u_dof)); + const auto K_up = + linear_operator(tangent_matrix.block(u_dof, p_dof)); + const auto K_pu = + linear_operator(tangent_matrix.block(p_dof, u_dof)); + const auto K_Jp = + linear_operator(tangent_matrix.block(J_dof, p_dof)); + const auto K_JJ = + linear_operator(tangent_matrix.block(J_dof, J_dof)); + PreconditionSelector, Vector> + preconditioner_K_Jp_inv("jacobi"); + preconditioner_K_Jp_inv.use_matrix( + tangent_matrix.block(J_dof, p_dof)); + ReductionControl solver_control_K_Jp_inv( + tangent_matrix.block(J_dof, p_dof).m() * + parameters.max_iterations_lin, + 1.0e-30, + parameters.tol_lin); + SolverSelector> solver_K_Jp_inv; + solver_K_Jp_inv.select("cg"); + solver_K_Jp_inv.set_control(solver_control_K_Jp_inv); + const auto K_Jp_inv = + inverse_operator(K_Jp, solver_K_Jp_inv, preconditioner_K_Jp_inv); + const auto K_pJ_inv = transpose_operator(K_Jp_inv); + const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv; + const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu; + const auto K_uu_con = K_uu + K_uu_bar_bar; + PreconditionSelector, Vector> + preconditioner_K_con_inv(parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner_K_con_inv.use_matrix( + tangent_matrix.block(u_dof, u_dof)); + ReductionControl solver_control_K_con_inv( + tangent_matrix.block(u_dof, u_dof).m() * + parameters.max_iterations_lin, + 1.0e-30, + parameters.tol_lin); + SolverSelector> solver_K_con_inv; + solver_K_con_inv.select("cg"); + solver_K_con_inv.set_control(solver_control_K_con_inv); + const auto K_uu_con_inv = + inverse_operator(K_uu_con, + solver_K_con_inv, + preconditioner_K_con_inv); + d_u = + K_uu_con_inv * (f_u - K_up * (K_Jp_inv * f_J - K_pp_bar * f_p)); + timer.leave_subsection(); + timer.enter_subsection("Linear solver postprocessing"); + std::cout << " PP " << std::flush; + d_J = K_pJ_inv * (f_p - K_pu * d_u); + d_p = K_Jp_inv * (f_J - K_JJ * d_J); + lin_it = solver_control_K_con_inv.last_step(); + lin_res = solver_control_K_con_inv.last_value(); + } + else if (parameters.type_lin == "Direct") + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(tangent_matrix); + A_direct.vmult(newton_update, system_rhs); + lin_it = 1; + lin_res = 0.0; + std::cout << " -- " << std::flush; + } + else + Assert(false, ExcMessage("Linear solver type not implemented")); + timer.leave_subsection(); + constraints.distribute(newton_update); + } + return std::make_pair(lin_it, lin_res); + } + template + void + Solid::output_results() const + { + DataOut data_out; + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + std::vector solution_name(dim, "displacement"); + solution_name.push_back("pressure"); + solution_name.push_back("dilatation"); + data_out.attach_dof_handler(dof_handler_ref); + data_out.add_data_vector(solution_n, + solution_name, + DataOut::type_dof_data, + data_component_interpretation); + Vector soln(solution_n.size()); + for (unsigned int i = 0; i < soln.size(); ++i) + soln(i) = solution_n(i); + MappingQEulerian q_mapping(degree, dof_handler_ref, soln); + data_out.build_patches(q_mapping, degree); + std::ostringstream filename; + filename << "solution-" << dim << "d-" << time.get_timestep() << ".vtk"; + std::ofstream output(filename.str().c_str()); + data_out.write_vtk(output); + } +} // namespace Step44 +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + using namespace Step44; + try + { + const unsigned int dim = 3; + Solid solid(SOURCE_DIR "/prm/parameters-step-44.prm"); + solid.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +} diff --git a/tests/adolc/step-44.output b/tests/adolc/step-44.output new file mode 100644 index 0000000000..9f2261a1cb --- /dev/null +++ b/tests/adolc/step-44.output @@ -0,0 +1,11 @@ + +DEAL::Timestep 1: 0.00 -0.000129 0.00 +DEAL::Timestep 2: 0.00 -0.000256 0.00 +DEAL::Timestep 3: 0.00 -0.000376 0.00 +DEAL::Timestep 4: 0.00 -0.000485 0.00 +DEAL::Timestep 5: 0.00 -0.000580 0.00 +DEAL::Timestep 6: 0.00 -0.000660 0.00 +DEAL::Timestep 7: 0.00 -0.000727 0.00 +DEAL::Timestep 8: 0.00 -0.000783 0.00 +DEAL::Timestep 9: 0.00 -0.000828 0.00 +DEAL::Timestep 10: 0.00 -0.000866 0.00 diff --git a/tests/adolc/symmetric_tensor_functions_02_1.cc b/tests/adolc/symmetric_tensor_functions_02_1.cc new file mode 100644 index 0000000000..c3ba107d2c --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_02_1.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Tensor inverse +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/symmetric_tensor_functions_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/symmetric_tensor_functions_02_1.output b/tests/adolc/symmetric_tensor_functions_02_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_02_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/symmetric_tensor_functions_02_2.cc b/tests/adolc/symmetric_tensor_functions_02_2.cc new file mode 100644 index 0000000000..259532bd02 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_02_2.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Tensor inverse +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/symmetric_tensor_functions_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/symmetric_tensor_functions_02_2.output b/tests/adolc/symmetric_tensor_functions_02_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_02_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/symmetric_tensor_functions_03_1.cc b/tests/adolc/symmetric_tensor_functions_03_1.cc new file mode 100644 index 0000000000..03a52ad7ab --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_03_1.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Eigenvalues +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/symmetric_tensor_functions_03.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output b/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output b/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output @@ -0,0 +1 @@ + diff --git a/tests/adolc/symmetric_tensor_functions_03_2.cc b/tests/adolc/symmetric_tensor_functions_03_2.cc new file mode 100644 index 0000000000..a0a74dd7ff --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_03_2.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Eigenvalues +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/symmetric_tensor_functions_03.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/symmetric_tensor_functions_03_2.output b/tests/adolc/symmetric_tensor_functions_03_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_03_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/symmetric_tensor_functions_04_1.cc b/tests/adolc/symmetric_tensor_functions_04_1.cc new file mode 100644 index 0000000000..735cf5811a --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_04_1.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Eigenvalues and eignvectors +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/symmetric_tensor_functions_04.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output b/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output b/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output @@ -0,0 +1 @@ + diff --git a/tests/adolc/symmetric_tensor_functions_04_2.cc b/tests/adolc/symmetric_tensor_functions_04_2.cc new file mode 100644 index 0000000000..1b03ca8751 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_04_2.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Eigenvalues and eignvectors +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/symmetric_tensor_functions_04.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/symmetric_tensor_functions_04_2.output b/tests/adolc/symmetric_tensor_functions_04_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/symmetric_tensor_functions_04_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/tensor_functions_02_1.cc b/tests/adolc/tensor_functions_02_1.cc new file mode 100644 index 0000000000..019aca9eb7 --- /dev/null +++ b/tests/adolc/tensor_functions_02_1.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Tensor inverse +// +// AD number type: ADOL-C taped + +#include "../ad_common_tests/tensor_functions_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor<2, double, AD::NumberTypes::adolc_taped>(); + test_tensor<3, double, AD::NumberTypes::adolc_taped>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/tensor_functions_02_1.output b/tests/adolc/tensor_functions_02_1.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/tensor_functions_02_1.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/tensor_functions_02_2.cc b/tests/adolc/tensor_functions_02_2.cc new file mode 100644 index 0000000000..632453f40c --- /dev/null +++ b/tests/adolc/tensor_functions_02_2.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test to check that tensor functions both compile and produce the right +// result when differentiated using the various auto-differentiable number +// types: Tensor inverse +// +// AD number type: ADOL-C tapeless + +#include "../ad_common_tests/tensor_functions_02.h" +#include "../tests.h" + +int +main() +{ + initlog(); + + deallog.push("Double"); + { + test_tensor<2, double, AD::NumberTypes::adolc_tapeless>(); + test_tensor<3, double, AD::NumberTypes::adolc_tapeless>(); + deallog << "OK" << std::endl; + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/tensor_functions_02_2.output b/tests/adolc/tensor_functions_02_2.output new file mode 100644 index 0000000000..3c6dd257b1 --- /dev/null +++ b/tests/adolc/tensor_functions_02_2.output @@ -0,0 +1,3 @@ + +DEAL::Double:OK +DEAL::OK diff --git a/tests/adolc/tensor_ops_01.cc b/tests/adolc/tensor_ops_01.cc new file mode 100644 index 0000000000..11295a53fd --- /dev/null +++ b/tests/adolc/tensor_ops_01.cc @@ -0,0 +1,236 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// compilation check constructor and operations acting on tensors of ADOL-C +// number types + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test_tensor() +{ + typedef typename Differentiation::AD::NumberTraits::ad_type + ad_number_t; + typedef Tensor<2, dim, ad_number_t> AD_Tensor; + typedef Tensor<2, dim, number_t> NonAD_Tensor; + + // Constructors + NonAD_Tensor t1; + AD_Tensor adt1, adt2; + AD_Tensor adt3(adt1); + + // Assignment + for (unsigned int i = 0; i < dim; ++i) + { + adt1[i][i] += 1.0; + for (unsigned int j = 0; j < dim; ++j) + { + t1[i][j] = dim * i + j; + adt1[i][j] += 2 * dim * i + j + 1; + adt2[i][j] += 3 * dim * i + j + 2; + } + } + AD_Tensor adt4 = adt1; + + // Comparison + const bool res_1 = (adt1 == adt2); + const bool res_2 = (adt1 != adt2); + + // Scalar operations + adt4 *= 2.0; + adt4 /= 4.0; + adt4[0][0] += 1.0; + adt4[0][0] *= 3.0; + + // Vector space operations + adt4 = 2.0 * adt1; + adt4 += adt2; + adt4 -= adt2; + adt4 += t1; + adt4 -= t1; + const AD_Tensor adt5 = adt1 + adt2; + const AD_Tensor adt6 = adt1 - adt2; + + // Contractions + const ad_number_t ad_res1 = scalar_product(adt1, adt2); + const ad_number_t ad_res2 = double_contract(adt1, adt2); + // ad_res1 = double_contract(adt1,t1); // TODO: Not defined. Conflicting + // number types ad_res1 = double_contract(t1,adt2); // TODO: Not defined. + // Conflicting number types + AD_Tensor adt7 = adt1 * adt2; + adt7 = t1 * adt2; + adt7 = adt1 * t1; + + // Outer product + const Tensor<1, dim, number_t> v1; + const Tensor<1, dim, ad_number_t> av1; + const AD_Tensor adt8 = outer_product(v1, av1); + const AD_Tensor adt9 = outer_product(av1, v1); + const AD_Tensor adt10 = outer_product(av1, av1); + + // Special operations + const ad_number_t det = determinant(adt1); + const ad_number_t tr = trace(adt1); + // const ad_number_t l1 = l1_norm(adt1); // TODO: Defined with std:: + // math operator const ad_number_t linf = linfty_norm(adt1); // TODO: Defined + // with std:: math operator + const AD_Tensor inv = invert(adt1); + const AD_Tensor trans = transpose(adt1); +} + +template +void +test_symmetric_tensor() +{ + typedef typename Differentiation::AD::NumberTraits::ad_type + ad_number_t; + typedef SymmetricTensor<2, dim, ad_number_t> AD_STensor; + typedef SymmetricTensor<4, dim, ad_number_t> AD_STensor4; + typedef SymmetricTensor<2, dim, number_t> NonAD_STensor; + typedef SymmetricTensor<4, dim, number_t> NonAD_STensor4; + + // Constructors + NonAD_STensor t1; + NonAD_STensor4 t2; + AD_STensor adt1, adt2; + AD_STensor adt3(adt1); + + // Assignment + for (unsigned int i = 0; i < dim; ++i) + { + for (unsigned int j = i; j < dim; ++j) + { + t1[i][j] = dim * i + j; + adt1[i][j] = 2 * dim * i + j + 1; + adt2[i][j] += 3 * dim * i + j + 2; + } + } + AD_STensor adt4 = adt1; + + // Comparison + bool res_1 = (adt1 == adt2); + bool res_2 = (adt1 != adt2); + + // Scalar operations + adt4 *= 2.0; + adt4 /= 4.0; + adt4 = 2.0 * adt1; + adt4[0][0] += 1.0; + adt4[0][0] -= 0.5; + adt4[0][0] *= 3.0; + adt4[0][0] /= 1.5; + + // Vector space operations + adt4 += adt2; + adt4 -= adt2; + adt4 += t1; + adt4 -= t1; + const AD_STensor adt5 = adt1 + adt2; + const AD_STensor adt6 = adt1 - adt2; + + // Contractions + const ad_number_t ad_res1 = scalar_product(adt1, adt2); + const ad_number_t ad_res2 = adt1 * adt2; + // const ad_number_t ad_res3 = scalar_product(t1,adt2); + // const ad_number_t ad_res4 = scalar_product(adt1,t1); + // const ad_number_t ad_res5 = t1*adt2; // TODO: Overload issue + // const ad_number_t ad_res6 = adt1*t1; // TODO: Overload issue + + // Outer product + const Tensor<1, dim, number_t> v1; + const Tensor<1, dim, ad_number_t> av1; + const AD_STensor adt8 = symmetrize(outer_product(v1, av1)); + const AD_STensor adt9 = symmetrize(outer_product(av1, v1)); + const AD_STensor adt10 = symmetrize(outer_product(av1, av1)); + + // Tensor outer product + const AD_STensor4 adt11 = outer_product(adt1, adt2); + // const AD_STensor4 adt12 = outer_product(t1,adt1); // TODO: Cannot deduce + // return type + + // Special operations + const ad_number_t det = determinant(adt1); + const ad_number_t tr = trace(adt1); + const ad_number_t inv1 = first_invariant(adt1); + const ad_number_t inv2 = second_invariant(adt1); + const ad_number_t inv3 = third_invariant(adt1); // TODO: Return type incorrect + // const ad_number_t l1 = l1_norm(adt1); // TODO: Defined with std:: math + // operator const ad_number_t linf = linfty_norm(adt1); // TODO: Defined with + // std:: math operator + const AD_STensor inv = invert(adt1); + const AD_STensor trans = transpose(adt1); + const AD_STensor dev = deviator(adt1); + + // Special tensors + const SymmetricTensor<2, dim, ad_number_t> ust = + unit_symmetric_tensor(); + const SymmetricTensor<4, dim, ad_number_t> dt = + deviator_tensor(); + // SymmetricTensor<4,dim,ad_number_t> idt = + // identity_tensor<4,dim,ad_number_t> (); // TODO: No number type + const SymmetricTensor<4, dim, ad_number_t> op = outer_product(adt1, adt2); + // SymmetricTensor<4,dim,ad_number_t> inv_op = invert(op); // TODO: Not + // implemented! +} + +int +main() +{ + initlog(); + + // It is assumed that whatever works for rank-2 + // tensors will work for the other ranks + + // --- Taped --- + // double type + test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor<3, + double, + Differentiation::AD::NumberTypes::adolc_taped>(); + // complex double + // test_tensor<3,std::complex,Differentiation::AD::NumberTypes::adolc_taped>(); + // test_symmetric_tensor<3,std::complex,Differentiation::AD::NumberTypes::adolc_taped>(); + + // --- Tapeless --- + // double type + test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor<3, + double, + Differentiation::AD::NumberTypes::adolc_tapeless>(); + // complex double + // test_tensor<3,std::complex,Differentiation::AD::NumberTypes::adolc_tapeless>(); + // test_symmetric_tensor<3,std::complex,Differentiation::AD::NumberTypes::adolc_tapeless>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/adolc/tensor_ops_02.cc b/tests/adolc/tensor_ops_02.cc new file mode 100644 index 0000000000..020da8ad5d --- /dev/null +++ b/tests/adolc/tensor_ops_02.cc @@ -0,0 +1,105 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check ability to cast from a tensor of ADOL-C numbers to floats +// In particular, this checks that the specializations for +// internal::NumberType::value() work as expected when +// Number == double and the input to value() is an ADOL-C number. + +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test_tensor() +{ + typedef typename Differentiation::AD::NumberTraits::ad_type + ad_number_t; + typedef Tensor<2, dim, ad_number_t> AD_Tensor; + typedef Tensor<2, dim, number_t> NonAD_Tensor; + + AD_Tensor adt1; + for (unsigned int i = 0; i < NonAD_Tensor::n_independent_components; ++i) + adt1[NonAD_Tensor::unrolled_to_component_indices(i)] = i + 1; + + const NonAD_Tensor t1 = NonAD_Tensor(adt1); + + Assert(t1.norm() > 0.0, ExcMessage("Cast and copy unsuccessful")); +} + +template +void +test_symmetric_tensor() +{ + typedef typename Differentiation::AD::NumberTraits::ad_type + ad_number_t; + typedef SymmetricTensor<2, dim, ad_number_t> AD_STensor2; + typedef SymmetricTensor<4, dim, ad_number_t> AD_STensor4; + typedef SymmetricTensor<2, dim, number_t> NonAD_STensor2; + typedef SymmetricTensor<4, dim, number_t> NonAD_STensor4; + + // Constructors + AD_STensor2 adt1; + AD_STensor4 adt2; + + for (unsigned int i = 0; i < NonAD_STensor2::n_independent_components; ++i) + adt1[NonAD_STensor2::unrolled_to_component_indices(i)] = i + 1; + + for (unsigned int i = 0; i < NonAD_STensor4::n_independent_components; ++i) + adt2.access_raw_entry(i) = i + 1; + + const NonAD_STensor2 t1 = NonAD_STensor2(adt1); + const NonAD_STensor4 t2 = NonAD_STensor4(adt2); + + Assert(t1.norm() > 0.0, ExcMessage("Cast and copy unsuccessful")); + Assert(t2.norm() > 0.0, ExcMessage("Cast and copy unsuccessful")); +} + +int +main() +{ + initlog(); + + // It is assumed that whatever works for rank-2 + // tensors will work for the other ranks + + // --- Taped --- + test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_taped>(); + test_symmetric_tensor<3, + double, + Differentiation::AD::NumberTypes::adolc_taped>(); + + // --- Tapeless --- + test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_tapeless>(); + test_symmetric_tensor<3, + double, + Differentiation::AD::NumberTypes::adolc_tapeless>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/adolc/tensor_ops_02.output b/tests/adolc/tensor_ops_02.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/adolc/tensor_ops_02.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5