Error using ADCZM

[chunhuizhao478]

In the computeInterfaceTraction() function, I try to approximate rate using current and old interface_displacement_jump as follows:

//Approximate slip rate
  RealVectorValue interface_displacement_jump_rate = (1.0/_dt) * (_interface_displacement_jump[_qp] - _interface_displacement_jump_old[_qp]) + Vini_vec;

However, I got the following error:

/Users/andyz/mambaforge3/envs/moose/libmesh/include/libmesh/type_vector.h:548:18: error: assigning to 'double' from incompatible type 'const MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>'
    _coords[i] = p._coords[i];
                 ^~~~~~~~~~~~
/Users/andyz/mambaforge3/envs/moose/libmesh/include/libmesh/vector_value.h:209:3: note: in instantiation of function template specialization 'libMesh::TypeVector<double>::TypeVector<MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>>' requested here
  TypeVector<T> (p)
  ^
/Users/andyz/projects/farms_rsf/src/materials/2d/implicit/AD/ADImplicitRateStateFrictionLaw2D.C:59:55: note: in instantiation of function template specialization 'libMesh::VectorValue<double>::VectorValue<MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>>' requested here
  RealVectorValue interface_displacement_jump_rate =  (_interface_displacement_jump[_qp] - _interface_displacement_jump_old[_qp]) + Vini_vec;
                                                      ^
1 error generated.
make: *** [/Users/andyz/projects/moose/framework/build.mk:149: /Users/andyz/projects/farms_rsf/build/unity_src/materials_2d_implicit_AD_Unity.arm64-apple-darwin20.0.0.opt.lo] Error 1

Any idea on how to get rid of this? Thanks a lot

[GiudGiud]

_coords[i] = MetaPhysicL::raw_value(p._coords[i]);

or

_coords[i] = p._coords[i].value();

this is happening because you are copying an AD number into a regular floating point number, hence removing the derivatives.
We do want to warn devs about removing derivatives like this

[GiudGiud]

no this is allowed. You just need to create the non-AD vector before computing the operations on the vectors and setting the non-AD vector

[GiudGiud]

you may be able to use raw_value on _interface_displacement_jump[_qp] directly actually

[chunhuizhao478]

sorry, could you write the format down? I'm not familiar with this AD, thanks!

[GiudGiud]

does this work?
either:

  RealVectorValue interface_displacement_jump_rate = (1.0/_dt) * (MetaPhysicL::raw_value(_interface_displacement_jump[_qp]) - _interface_displacement_jump_old[_qp]) + Vini_vec;

or using AD (probably preferrable actually)

  ADRealVectorValue interface_displacement_jump_rate = (1.0/_dt) * (_interface_displacement_jump[_qp] - _interface_displacement_jump_old[_qp]) + Vini_vec;

[chunhuizhao478]

Thanks! The first one works.

[GiudGiud]

ok first one loses the derivatives so just make sure you dont use it to compute an ADReal residual!

[chunhuizhao478]

Eventually the value will be put into _interface_traction(declareADPropertyByName<RealVectorValue>(_base_name + "interface_traction")), which is a RealVectorValue, I think it's ok?

[GiudGiud]

Nope not ok.
It’s an AD material property, the type is actually ADRealVectorValue

[chunhuizhao478]

The error message is as follows (I have updated moose to latest version):

(moose) andyz@Chunhuis-MBP-2 farms_rsf % make -j8
Compiling C++ (in opt mode) /Users/andyz/projects/farms_rsf/build/unity_src/materials_2d_implicit_Unity.C...
In file included from /Users/andyz/projects/farms_rsf/build/unity_src/materials_2d_implicit_Unity.C:2:
/Users/andyz/projects/farms_rsf/src/materials/2d/implicit/ADImplicitRateStateFrictionLaw2D.C:71:23: error: no matching function for call to 'asinh'
  ADReal mu = rsf_a * asinh( sliprate_strike/(2*delta_o) * exp((f_o + rsf_b * log(delta_o * statevar/rsf_L))/rsf_a) );
                      ^~~~~
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/math.h:341:15: note: candidate function not viable: no known conversion from 'typename MultipliesType<DualNumber<double, SemiDynamicSparseNumberArray<double, unsigned long long, NWrapper<53>>, true>, double, false>::supertype' (aka 'DualNumber<double, SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>') to 'double' for 1st argument
extern double asinh(double);
              ^
/Users/andyz/mambaforge3/envs/moose/bin/../include/c++/v1/math.h:1093:46: note: candidate function not viable: no known conversion from 'typename MultipliesType<DualNumber<double, SemiDynamicSparseNumberArray<double, unsigned long long, NWrapper<53>>, true>, double, false>::supertype' (aka 'DualNumber<double, SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>') to 'float' for 1st argument
inline _LIBCPP_INLINE_VISIBILITY float       asinh(float __lcpp_x) _NOEXCEPT       {return ::asinhf(__lcpp_x);}
                                             ^
/Users/andyz/mambaforge3/envs/moose/bin/../include/c++/v1/math.h:1094:46: note: candidate function not viable: no known conversion from 'typename MultipliesType<DualNumber<double, SemiDynamicSparseNumberArray<double, unsigned long long, NWrapper<53>>, true>, double, false>::supertype' (aka 'DualNumber<double, SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>') to 'long double' for 1st argument
inline _LIBCPP_INLINE_VISIBILITY long double asinh(long double __lcpp_x) _NOEXCEPT {return ::asinhl(__lcpp_x);}
                                             ^
/Users/andyz/mambaforge3/envs/moose/bin/../include/c++/v1/math.h:1099:1: note: candidate template ignored: requirement 'std::is_integral<MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>>::value' was not satisfied [with _A1 = MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long long, MetaPhysicL::NWrapper<53>>, true>]
asinh(_A1 __lcpp_x) _NOEXCEPT {return ::asinh((double)__lcpp_x);}
^
1 error generated.
make: *** [/Users/andyz/projects/moose/framework/build.mk:145: /Users/andyz/projects/farms_rsf/build/unity_src/materials_2d_implicit_Unity.arm64-apple-darwin20.0.0.opt.lo] Error 1
(moose) andyz@Chunhuis-MBP-2 farms_rsf % 

the header file:

#pragma once

#include "ADCZMComputeLocalTractionTotalBase.h"

class ADImplicitRateStateFrictionLaw2D : public ADCZMComputeLocalTractionTotalBase
{
public:
  static InputParameters validParams();
  ADImplicitRateStateFrictionLaw2D(const InputParameters & parameters);

protected:
  virtual void computeInterfaceTraction();

  //rate-and-state friction coefficients
  Real _f_o;
  Real _rsf_a;
  Real _rsf_b;
  Real _rsf_L;
  Real _delta_o;
  Real _Vini;

  //initial shear and normal traction
  Real _Tn_o;
  Real _Ts_o;

  //old/older interface jump
  const MaterialProperty<RealVectorValue> & _interface_displacement_jump_old;
  const MaterialProperty<RealVectorValue> & _interface_displacement_jump_older;
};

the source file:

#include "ADImplicitRateStateFrictionLaw2D.h"

registerMooseObject("farms_rsfApp", ADImplicitRateStateFrictionLaw2D);

InputParameters
ADImplicitRateStateFrictionLaw2D::validParams()
{
  InputParameters params = ADImplicitRateStateFrictionLaw2D::validParams();
  return params;
}

ADImplicitRateStateFrictionLaw2D::ADImplicitRateStateFrictionLaw2D(const InputParameters & parameters)
  : ADCZMComputeLocalTractionTotalBase(parameters),
  _f_o(getParam<Real>("f_o")),
  _rsf_a(getParam<Real>("rsf_a")),
  _rsf_b(getParam<Real>("rsf_b")),
  _rsf_L(getParam<Real>("rsf_L")),
  _delta_o(getParam<Real>("delta_o")),
  _Vini(getParam<Real>("Vini")),
  _Tn_o(getParam<Real>("Tn_o")),
  _Ts_o(getParam<Real>("Ts_o")),
  _interface_displacement_jump_old(getMaterialPropertyOldByName<RealVectorValue>(_base_name + "interface_displacement_jump")),
  _interface_displacement_jump_older(getMaterialPropertyOlderByName<RealVectorValue>(_base_name + "interface_displacement_jump"))
{
}

//Define Frictional Law (function of sliprate and statevar)
double mu_friction_law_implicit2DAD(Real sliprate, Real statevar, Real rsf_a, Real rsf_b, Real rsf_L, Real delta_o, Real f_o)
{
  double mu = 0;
  mu = rsf_a * asinh( sliprate/(2*delta_o) * exp((f_o + rsf_b * log(delta_o * statevar/rsf_L))/rsf_a) );
  return mu;
}

void
ADImplicitRateStateFrictionLaw2D::computeInterfaceTraction()
{
  //Define Parameters
  Real f_o = _f_o;
  Real rsf_a = _rsf_a;
  Real rsf_b = _rsf_b;
  Real rsf_L = _rsf_L;
  Real delta_o = _delta_o;
  Real Vini = _Vini;
  Real Tn_o = _Tn_o;
  Real Ts_o = _Ts_o;

  //Get initial slip rate value
  RealVectorValue Vini_vec(0.0,Vini,0.0);

  //Get normal traction
  //if (_interface_displacement_jump[_qp](0) < 0){}else{}
  Real Tn = abs(Tn_o);

  //Get old state variable
  Real statevar_old = 0;

  //Approximate slip rate
  ADRealVectorValue interface_displacement_jump_rate = (1.0/_dt) * (_interface_displacement_jump[_qp] - _interface_displacement_jump_old[_qp]) + Vini_vec;
  RealVectorValue interface_displacement_jump_rate_old = (1.0/_dt) * (_interface_displacement_jump_old[_qp] - _interface_displacement_jump_older[_qp]) + Vini_vec;

  //Get strike direction at t
  ADReal sliprate_strike = interface_displacement_jump_rate(1);
  Real sliprate_strike_old = interface_displacement_jump_rate_old(1);

  //Define state variable at t+dt
  Real statevar = (rsf_L / sliprate_strike_old) + (statevar_old - rsf_L / sliprate_strike_old) * exp( -1 * sliprate_strike_old * _dt / rsf_L );
  
  //Define friction mu at t+dt
  //!arcsinh(z)= ln(z+sqrt(z^2+1))
  ADReal mu = rsf_a * asinh( sliprate_strike/(2*delta_o) * exp((f_o + rsf_b * log(delta_o * statevar/rsf_L))/rsf_a) );
  
  //Compute shear traction at t+dt
  ADReal Ts = Tn * mu;

  //Save results
  _statevar[_qp] = statevar;
  _sliprate[_qp] = sliprate_strike;

  //Feed traction into the system
  //Assign back traction in CZM
  ADRealVectorValue traction;

  traction(0) = 0.0; 
  traction(1) = -Ts+Ts_o; //"-"
  traction(2) = 0.0;

  _interface_traction[_qp] = traction;
}

[roystgnr]

I'm a little perturbed that the error message suggests it's not even considering the definition from dualnumber.h. Is there any possibility that the headers here ended up including dualnumber_forward.h but not dualnumber.h (or at least dualnumber_decl.h)?

[ahzeeshan]

I got the same error when using the log function from <cmath> but std::log works without errors in an ADKernel.

[lindsayad]

@roystgnr do you think if we put our dual number math function overloads in the MetaPhysicL namespace instead of std that argument dependent lookup would work for us here? Unfortunately doing that migration would break a lot of existing code which uses std::

[roystgnr]

It wouldn't be crazy to put them in both places. Backwards compatibility and (I hope; I'm not sure) ADL too.

[roystgnr]

Though I think using namespace std or explicit std::asinh etc. has to be the right thing to do here. Overloads for use with custom classes is supposed to be the use case where it's expected to put user code in namespace std.