Help Developing a Stress Induced Diffusion Kernel

[Supernova963]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hello,

I am trying to implement the following equation for stress induced diffusion,

The terms affected by the divergence (the dot product is missing), are the ones I am interested in implementing. As of now I have chosen to implement the first term via the matDiffusion kernel and am trying to create the second term in what I am calling the StressDrift kernel. I currently have the files made, but am unsure if the equations are implemented correctly as it is my first attempt at creating a new kernel.

Non-linear variable: $c^o$
Prefactor terms: $D, R, T, \eta, V$ (Input as a derivative material)
And I want to take the gradient of hydrostatic stress. (First order Lagrange)

I have directly based the code off matDiffusion as it seemed the most similar to what I want to do.

If anyone has any advice or feedback on whether it looks like it implements the equation please let me know. Anything is greatly appreciated.

StressDriftBase.h

#pragma once

#include "GenericKernelGrad.h"
#include "JvarMapInterface.h"
#include "DerivativeMaterialInterface.h"
#include "Coupleable.h"

template <typename T, bool is_ad>
class StressDriftBaseTempl
  : public DerivativeMaterialInterface<JvarMapKernelInterface<GenericKernelGrad<is_ad>>>
{
public:
  static InputParameters validParams();

  StressDriftBaseTempl(const InputParameters & parameters);

protected:
  virtual GenericRealVectorValue<is_ad> precomputeQpResidual() override;

  /// Prefactor
  const GenericMaterialProperty<T, is_ad> & _prefactor;

  /// Gradient of the stress
  const GenericVariableGradient<is_ad> & _grad_stress;

  usingGenericKernelGradMembers;
};

template <typename T>
class StressDriftBase : public StressDriftBaseTempl<T, false>
{
public:
  static InputParameters validParams();

  StressDriftBase(const InputParameters & parameters);

  using StressDriftBaseTempl<T, false>::StressDriftBaseTempl;

protected:
  virtual void initialSetup() override;
  virtual RealGradient precomputeQpJacobian() override;
  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
  virtual RealGradient precomputeQpCJacobian();

  /// prefector derivative w.r.t. the kernel variable
  const MaterialProperty<T> & _dprefactor_dc;

  /// prefactor derivatives w.r.t. coupled variables
  std::vector<const MaterialProperty<T> *> _dprefactor_darg;

  /// is the kernel used in a coupled form?
  const bool _is_coupled;

  /// int label for the Concentration
  unsigned int _v_var;

  /// int label for the stress
  unsigned int _stress_var;

  using StressDriftBaseTempl<T, false>::_prefactor;
  using StressDriftBaseTempl<T, false>::_qp;
  using StressDriftBaseTempl<T, false>::_grad_phi;
  using StressDriftBaseTempl<T, false>::_j;
  using StressDriftBaseTempl<T, false>::_i;
  using StressDriftBaseTempl<T, false>::_phi;
  using StressDriftBaseTempl<T, false>::_grad_stress; // new
  using StressDriftBaseTempl<T, false>::_grad_test;
  using StressDriftBaseTempl<T, false>::_var;
  using StressDriftBaseTempl<T, false>::_coupled_moose_vars;
};

StressDriftBase.C

#include "StressDriftBase.h"

template <typename T, bool is_ad>
InputParameters
StressDriftBaseTempl<T, is_ad>::validParams()
{
  InputParameters params = GenericKernelGrad<is_ad>::validParams();
  params.addParam<MaterialPropertyName>(
      "prefactor", "prefactor", "The prefactor value or material property");
  params.addCoupledVar("v",
                       "Coupled concentration variable for kernel to operate on; if this "
                       "is not specified, the kernel's nonlinear variable will be used as "
                       "usual");
  params.addCoupledVar("stress", "Coupled stress variable.");
  return params;
}

template <typename T>
InputParameters
StressDriftBase<T>::validParams()
{
  InputParameters params = StressDriftBaseTempl<T, false>::validParams();
  params.addCoupledVar("args",
                       "Optional vector of arguments for the prefactor. If provided and "
                       "prefactor is a derivative parsed material, Jacobian contributions from "
                       "the prefactor will be automatically computed");
  return params;
}

template <typename T, bool is_ad>
StressDriftBaseTempl<T, is_ad>::StressDriftBaseTempl(const InputParameters & parameters)
  : DerivativeMaterialInterface<JvarMapKernelInterface<GenericKernelGrad<is_ad>>>(parameters),
    _prefactor(this->template getGenericMaterialProperty<T, is_ad>("prefactor")),
    _grad_stress(this->template coupledGenericGradient<is_ad>("stress"))
{
}

template <typename T>
StressDriftBase<T>::StressDriftBase(const InputParameters & parameters)
  : StressDriftBaseTempl<T, false>(parameters),
    _dprefactor_dc(this->template getMaterialPropertyDerivative<T>("prefactor", _var.name())),
    _dprefactor_darg(_coupled_moose_vars.size()),
    _is_coupled(this->isCoupled("v")),
    _v_var(_is_coupled ? this->coupled("v") : _var.number()),
    _stress_var(_var.number())
{
  // fetch derivatives
  for (unsigned int i = 0; i < _dprefactor_darg.size(); ++i)
    _dprefactor_darg[i] = &this->template getMaterialPropertyDerivative<T>(
        "prefactor", _coupled_moose_vars[i]->name());
}

template <typename T, bool is_ad>
GenericRealVectorValue<is_ad>
StressDriftBaseTempl<T, is_ad>::precomputeQpResidual()
{
  return _prefactor[_qp] * _grad_stress[_qp] ;
}

template <typename T>
void
StressDriftBase<T>::initialSetup()
{
  this->template validateNonlinearCoupling<T>("prefactor");
}

template <typename T>
RealGradient
StressDriftBase<T>::precomputeQpCJacobian()
{
  return _prefactor[_qp] * _grad_phi[_j][_qp];
}

template <typename T>
RealGradient
StressDriftBase<T>::precomputeQpJacobian()
{
  auto sum = _dprefactor_dc[_qp] * _phi[_j][_qp] * _grad_stress[_qp];
  if (!_is_coupled)
    sum += precomputeQpCJacobian();

  return sum;
}

template <typename T>
Real
StressDriftBase<T>::computeQpOffDiagJacobian(unsigned int jvar)
{
  // get the coupled variable jvar is referring to
  const auto cvar = this->mapJvarToCvar(jvar);

  auto sum = (*_dprefactor_darg[cvar])[_qp] * _phi[_j][_qp] * _grad_stress[_qp];
  if (_v_var == jvar)
    sum += precomputeQpCJacobian();

  return sum * _grad_test[_i][_qp];
}

template class StressDriftBase<Real>;
template class StressDriftBase<RealTensorValue>;
template class StressDriftBaseTempl<Real, false>;
template class StressDriftBaseTempl<RealTensorValue, false>;
template class StressDriftBaseTempl<Real, true>;
template class StressDriftBaseTempl<RealTensorValue, true>;

[GiudGiud]

hello

so the second term, is it

div ( Dc0 / RT grad (tr(sigma))) ?

and then you integrated this by parts right?

I'd expect to see the prefactor * the variable (usually _u[_qp], with u set using coupledValue) * grad_stress

is the variable (c0) value baked into the prefactor?

[Supernova963]

Yes the prefactor has $c^o$ as a term which is the nonlinear variable. I didn't know how to handled it in the code. How it's written currently I have the derivativematerial having $c^o$ as a coupled variable.

But yes, I integrated $\frac{Dc^o}{RT}\nabla(tr(\sigma))$.