Issue with Custom kernel

[ashishdhole]

Hello,
I made this new kernel to incorporate the friction pressure as explained in
10.1088/0965-0393/20/6/065008
Ideally, after sometime the friction pressure should get equal to or greater than driving pressure and there should not be any grain boundary movement. but for some reason this is not happening. Is there any issue with the code? Can you have a look at it.
The kernel is a combination of two existing kernels ACGrGrPoly and ACInterface with an extra condition and a term I called friction_energy.

#include "AdvanceFrictionPressure.h"

registerMooseObject("pfubcApp", AdvanceFrictionPressure);

InputParameters
AdvanceFrictionPressure::validParams()
{
  InputParameters params = JvarMapKernelInterface<Kernel>::validParams();
  params.addClassDescription("Gradient energy Allen-Cahn Kernel");
  params.addParam<MaterialPropertyName>("mob_name", "L", "The mobility used with the kernel");
  params.addParam<MaterialPropertyName>("kappa_name", "kappa_op", "The kappa used with the kernel");
  params.addParam<bool>("variable_L",
                        true,
                        "The mobility is a function of any MOOSE variable (if false, L must be "
                        "constant over the entire domain!)");
  params.addCoupledVar("g", 0.0, "Friction pressure");
  params.addRequiredCoupledVarWithAutoBuild(
      "v", "var_name_base", "op_num", "Array of coupled variables");
  params.addParam<MaterialPropertyName>("sigma", "GB energy (sigma)");
  return params;
}

AdvanceFrictionPressure::AdvanceFrictionPressure(const InputParameters & parameters)
  : DerivativeMaterialInterface<JvarMapKernelInterface<Kernel>>(parameters),
    _L(getMaterialProperty<Real>("mob_name")),
    _kappa(getMaterialProperty<Real>("kappa_name")),
    _variable_L(getParam<bool>("variable_L")),
    _dLdop(getMaterialPropertyDerivative<Real>("mob_name", _var.name())),
    _d2Ldop2(getMaterialPropertyDerivative<Real>("mob_name", _var.name(), _var.name())),
    _dkappadop(getMaterialPropertyDerivative<Real>("kappa_name", _var.name())),
    _dLdarg(_n_args),
    _d2Ldargdop(_n_args),
    _d2Ldarg2(_n_args),
    _dkappadarg(_n_args),
    _gradarg(_n_args),
    _g(coupledValue("g")),
    _mu(getMaterialProperty<Real>("mu")),
    _op_num(coupledComponents("v")),
    _vals(coupledValues("v")),
    _vals_var(coupledIndices("v")),
    _sigma(getMaterialProperty<Real>("sigma")),
    _gamma(getMaterialProperty<Real>("gamma_asymm"))
{
  // Initialize derivatives and gradients for mobility, kappa, and coupled variables
  for (unsigned int i = 0; i < _n_args; ++i)
  {
    MooseVariable * ivar = _coupled_standard_moose_vars[i];
    const VariableName iname = ivar->name();

    if (iname == _var.name())
    {
      // Validate that the kernel variable is not specified in coupled variables
      if (isCoupled("args"))
        paramError("args",
                   "The kernel variable should not be specified in the coupled `args` parameter.");
      else
        paramError("coupled_variables",
                   "The kernel variable should not be specified in the coupled `coupled_variables` "
                   "parameter.");
    }

    // Assign material property derivatives and gradients
    _dLdarg[i] = &getMaterialPropertyDerivative<Real>("mob_name", i);
    _dkappadarg[i] = &getMaterialPropertyDerivative<Real>("kappa_name", i);
    _d2Ldargdop[i] = &getMaterialPropertyDerivative<Real>("mob_name", iname, _var.name());
    _gradarg[i] = &(ivar->gradSln());

    _d2Ldarg2[i].resize(_n_args);
    for (unsigned int j = 0; j < _n_args; ++j)
      _d2Ldarg2[i][j] = &getMaterialPropertyDerivative<Real>("mob_name", i, j);
  }
}

void
AdvanceFrictionPressure::initialSetup()
{
  // Validate coupling for material properties
  validateCoupling<Real>("mob_name");
  validateCoupling<Real>("kappa_name");
}

Real
AdvanceFrictionPressure::assignThisOp()
{
  return _u[_qp];
}

std::vector<Real>
AdvanceFrictionPressure::assignOtherOps()
{
  std::vector<Real> other_ops(_op_num);
  for (unsigned int s = 0; s < _op_num; ++s)
    other_ops[s] = (*_vals[s])[_qp];

  return other_ops;
}

RealGradient
AdvanceFrictionPressure::gradL()
{
  RealGradient gradient = _grad_u[_qp] * _dLdop[_qp];
  for (unsigned int i = 0; i < _n_args; ++i)
    gradient += (*_gradarg[i])[_qp] * (*_dLdarg[i])[_qp];
  return gradient;
}

RealGradient
AdvanceFrictionPressure::nablaLPsi()
{
  // Compute the gradient of L*psi using the product rule
  RealGradient gradient = _L[_qp] * _grad_test[_i][_qp];

  if (_variable_L)
    gradient += gradL() * _test[_i][_qp];

  return gradient;
}

RealGradient
AdvanceFrictionPressure::kappaNablaLPsi()
{
  return _kappa[_qp] * nablaLPsi();
}

Real
AdvanceFrictionPressure::computeQpResidual()
{
  Real op = assignThisOp();
  std::vector<Real> other_ops(_op_num);
  other_ops = assignOtherOps();
  Real g = _g[_qp];

  // Initialize sums for order parameters
  Real SumOPGj = 0.0;
  Real SumOPj = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumOPj += other_ops[i] * other_ops[i];
    SumOPGj += other_ops[i];
  }

  // Compute intermediate values
  Real inter = _grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp];
  Real op_cubed = op * op * op;
  Real f_loc = _mu[_qp] * (op_cubed - op + 2.0 * _gamma[_qp] * op * SumOPj);
  Real total = f_loc + inter;
  Real friction_energy = 3.0 * g * op * SumOPGj;

  // Compute sign of total
  const int sgn = (total > 0) - (total < 0);

  // Initialize residual result
  Real epsilon = 1e-5;
  Real Res = 0.0;
  Real total_abs = std::abs(total);
  Real friction_energy_abs = std::abs(friction_energy);
  Real max_abs = std::max(total_abs, friction_energy_abs);

  // Compute residual if condition is met
  if (total_abs - friction_energy_abs > epsilon * max_abs)
  {
    Real op_cubed = op * op * op;
    Real op_term = op_cubed - op + 2.0 * _gamma[_qp] * op * SumOPj;
    Real grad_u_term = _grad_u[_qp] * kappaNablaLPsi();
    Real sigma_term = _L[_qp] * sgn * (3.0 * g * op * SumOPGj);

    Res = (_L[_qp] * _mu[_qp] * op_term) + grad_u_term - sigma_term;
  }
  else
  {
    Res = 0;
  }

  return Res;
}

Real
AdvanceFrictionPressure::computeQpJacobian()
{
  // Assign primary and other order parameters
  Real op = assignThisOp();
  std::vector<Real> other_ops = assignOtherOps(); // Direct assignment
  const Real g = _g[_qp];

  // Initialize sums for order parameters
  Real SumOPGj = 0.0;
  Real SumOPj = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumOPj += other_ops[i] * other_ops[i];
    SumOPGj += other_ops[i];
  }

  // Compute intermediate values
  Real inter = _grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp];
  Real op_cubed = op * op * op;
  Real f_loc = _mu[_qp] * (op_cubed - op + 2.0 * _gamma[_qp] * op * SumOPj);
  Real total = f_loc + inter;
  Real friction_energy = 3.0 * g * op * SumOPGj;

  // Compute sign of total
  const int sgn = (total > 0) - (total < 0);

  // Compute the derivative of the gradient of the mobility
  RealGradient dsum =
      (_dkappadop[_qp] * _L[_qp] + _kappa[_qp] * _dLdop[_qp]) * _phi[_j][_qp] * _grad_test[_i][_qp];

  if (_variable_L)
  {
    RealGradient dgradL =
        _grad_phi[_j][_qp] * _dLdop[_qp] + _grad_u[_qp] * _phi[_j][_qp] * _d2Ldop2[_qp];

    for (unsigned int i = 0; i < _n_args; ++i)
      dgradL += (*_gradarg[i])[_qp] * _phi[_j][_qp] * (*_d2Ldargdop[i])[_qp];

    dsum += (_kappa[_qp] * dgradL + _dkappadop[_qp] * _phi[_j][_qp] * gradL()) * _test[_i][_qp];
  }

  Real epsilon = 1e-5;
  Real total_abs = std::abs(total);
  Real friction_energy_abs = std::abs(friction_energy);
  Real max_abs = std::max(total_abs, friction_energy_abs);

  Real Jac = 0.0; // Initial value

  if (total_abs - friction_energy_abs > epsilon * max_abs)
  {
    Real op_squared = op * op;
    Real phi_term = _phi[_j][_qp] * (3.0 * op_squared - 1.0 + 2.0 * _gamma[_qp] * SumOPj);
    Real grad_phi_term = _grad_phi[_j][_qp] * kappaNablaLPsi();
    Real grad_u_term = _grad_u[_qp] * dsum;
    Real sigma_term = _L[_qp] * sgn * _phi[_j][_qp] * (3.0 * g * SumOPGj);

    Jac = (_L[_qp] * _mu[_qp] * phi_term) + grad_phi_term + grad_u_term - sigma_term;
  }
  else
  {
    Jac = 0;
  }

  return Jac;
}

I deal behaviour

But I am getting something like this

Thank you

[GiudGiud]

Hello

so this term is supposed to be added only in certain conditions and counter-balance another term in the equations?

[GiudGiud]

Not really.
These are not very many qp were it is achieved either

[ashishdhole]

how to output the element number in this case, maybe that will help?
like _qp for quadrature point index, is there any defined variable for element number?

[GiudGiud]

_current_elem->id() gives you the element number

[ashishdhole]

So I can see that some elements do reach the condition

102450
102450
102451
102451
102451
102451
102451
102451
102451
102452
102452
102452
102454
102454
102454
102454
102454
102454
102460
102460
102460
102460
102460
102462
102462
102462
102462
102462
102462
 3 Nonlinear |R| = 5.185226e-04
 Solve Converged!
  Finished Solving                                                                       [123.66 s] [  531 MB]
  Currently Computing User Objects
    Finished Flooding Features                                                           [  5.10 s] [  531 MB]
  Finished Computing User Objects                                                        [  5.80 s] [  531 MB]

Outlier Variable Residual Norms:
  gr2: 2.556179e-04

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+
| time           | R              | avg_grain_vol  | g              | grain_tracker  |
+----------------+----------------+----------------+----------------+----------------+
|   0.000000e+00 |   0.000000e+00 |   1.883946e+02 |   0.000000e+00 |   2.000000e+02 |
|   5.000000e-02 |   7.821093e+00 |   1.921696e+02 |   3.120755e+05 |   2.000000e+02 |
+----------------+----------------+----------------+----------------+----------------+

Steady-State Relative Differential Norm: 1.21968

but each element is repeated 5 times. I don't know why that us.

I think I should play with the material properties to see if halt happens.
I will get back to you

[GiudGiud]

but each element is repeated 5 times. I don't know why that us.

probably 5 quadrature points. you could output both qp and the element id to check