Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Adolc tensor refactor #14205

Closed
wants to merge 6 commits into from
Closed

Adolc tensor refactor #14205

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
1675bc6
Fix abs() return types in l1_norm() and linfty_norm().
drwells Aug 6, 2022
db0906f
Combine both implementations of l1_norm() and linfty_norm().
drwells Aug 6, 2022
cf015cc
Actually test the l1_norm() and linfty_norm() with adol-c.
drwells Aug 6, 2022
1f85371
Remove unused headers.
drwells Aug 6, 2022
bbd2e45
Delete an unnecessary check for adol-c.
drwells Aug 7, 2022
efe392b
Clean up a test.
drwells Aug 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
79 changes: 10 additions & 69 deletions include/deal.II/base/tensor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,8 @@
#include <deal.II/base/exceptions.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/table_indices.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/tensor_accessors.h>

#ifdef DEAL_II_WITH_ADOLC
# include <adolc/adouble.h> // Taped double
#endif

#include <cmath>
#include <ostream>

Expand Down Expand Up @@ -1113,14 +1108,6 @@ template <typename OtherNumber>
constexpr inline bool
Tensor<0, dim, Number>::operator==(const Tensor<0, dim, OtherNumber> &p) const
{
# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
Assert(!(std::is_same<Number, adouble>::value ||
std::is_same<OtherNumber, adouble>::value),
ExcMessage(
"The Tensor equality operator for ADOL-C taped numbers has not yet "
"been extended to support advanced branching."));
# endif

return numbers::values_are_equal(value, p.value);
}

Expand Down Expand Up @@ -3003,17 +2990,18 @@ project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &A);
* @relatesalso Tensor
*/
template <int dim, typename Number>
inline Number
inline typename numbers::NumberTraits<Number>::real_type
l1_norm(const Tensor<2, dim, Number> &t)
{
Number max = internal::NumberType<Number>::value(0.0);
typename numbers::NumberTraits<Number>::real_type max(0.0);
for (unsigned int j = 0; j < dim; ++j)
{
Number sum = internal::NumberType<Number>::value(0.0);
typename numbers::NumberTraits<Number>::real_type sum(0.0);
for (unsigned int i = 0; i < dim; ++i)
sum += numbers::NumberTraits<Number>::abs(t[i][j]);

if (sum > max)
// For compatibility with ADOL-C we need the comparison written this way
if (sum - max > 0.0)
max = sum;
}

Expand All @@ -3029,17 +3017,18 @@ l1_norm(const Tensor<2, dim, Number> &t)
* @relatesalso Tensor
*/
template <int dim, typename Number>
inline Number
inline typename numbers::NumberTraits<Number>::real_type
linfty_norm(const Tensor<2, dim, Number> &t)
{
Number max = internal::NumberType<Number>::value(0.0);
typename numbers::NumberTraits<Number>::real_type max(0.0);
for (unsigned int i = 0; i < dim; ++i)
{
Number sum = internal::NumberType<Number>::value(0.0);
typename numbers::NumberTraits<Number>::real_type sum(0.0);
for (unsigned int j = 0; j < dim; ++j)
sum += numbers::NumberTraits<Number>::abs(t[i][j]);

if (sum > max)
// For compatibility with ADOL-C we need the comparison written this way
if (sum - max > 0.0)
max = sum;
}

Expand All @@ -3048,54 +3037,6 @@ linfty_norm(const Tensor<2, dim, Number> &t)

/** @} */


#ifndef DOXYGEN


# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING

// Specialization of functions for ADOL-C number types when
// the advanced branching feature is used
template <int dim>
inline adouble
l1_norm(const Tensor<2, dim, adouble> &t)
{
adouble max = internal::NumberType<adouble>::value(0.0);
for (unsigned int j = 0; j < dim; ++j)
{
adouble sum = internal::NumberType<adouble>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
sum += std::fabs(t[i][j]);

condassign(max, (sum > max), sum, max);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@drwells For these taped types, one needs to use the condassign() function to switch between branches. An if statement won't work if the evaluation point changes which branch is taken.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll write a test so that I understand this better.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had the ADOL-C manual in hand when decomposing and writing these functions, so whenever you come across an unusual implementation then I'd suggest consulting the manual. The requirements of the taping algorithm often dictates how code with branches has to be written in order to be taped AD compatible. It would not surprise me if there is not a test for every single specialised function -- the implementation was probably "best effort".

}

return max;
}


template <int dim>
inline adouble
linfty_norm(const Tensor<2, dim, adouble> &t)
{
adouble max = internal::NumberType<adouble>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
{
adouble sum = internal::NumberType<adouble>::value(0.0);
for (unsigned int j = 0; j < dim; ++j)
sum += std::fabs(t[i][j]);

condassign(max, (sum > max), sum, max);
}

return max;
}

# endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif
3 changes: 0 additions & 3 deletions include/deal.II/base/tensor_accessors.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,6 @@

#include <deal.II/base/config.h>

#include <deal.II/base/table_indices.h>
#include <deal.II/base/template_constraints.h>


DEAL_II_NAMESPACE_OPEN

Expand Down
20 changes: 10 additions & 10 deletions tests/adolc/tensor_ops_01.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ test_tensor()
// Comparison
const bool res_1 = (adt1 == adt2);
const bool res_2 = (adt1 != adt2);
const bool res_3 = adt1 == t1;
const bool res_4 = t1 == adt1;

// Scalar operations
adt4 *= 2.0;
Expand Down Expand Up @@ -97,13 +99,12 @@ test_tensor()
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);
const ad_number_t det = determinant(adt1);
const ad_number_t tr = trace(adt1);
const ad_number_t l1 = l1_norm(adt1);
const ad_number_t linf = linfty_norm(adt1);
const AD_Tensor inv = invert(adt1);
const AD_Tensor trans = transpose(adt1);
}

template <int dim,
Expand Down Expand Up @@ -185,9 +186,8 @@ test_symmetric_tensor()
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_number_t l1 = l1_norm(adt1);
// const ad_number_t linf = linfty_norm(adt1);
const AD_STensor inv = invert(adt1);
const AD_STensor trans = transpose(adt1);
const AD_STensor dev = deviator(adt1);
Expand Down
0 ...dolc_with_advanced_branching=false.output → tests/adolc/tensor_ops_01.output
View file
Open in desktop
File renamed without changes.
2 changes: 0 additions & 2 deletions tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output
View file
Open in desktop

This file was deleted.

2 changes: 0 additions & 2 deletions tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output
View file
Open in desktop

This file was deleted.