Difference in the pinning pressure results. MOOSE Vs. Analytical solutions

[ashishdhole]

Hello,

I am working on pinning pressure model following the paper (10.1088/0965-0393/20/6/065008)
where they have expressed the pinning pressure in the evolution equation as second term in the following equation

Now I am working on this kernel for some time and I finally have a kernel that debugger shows "no error!"
So I tried to compare with the analytical model
Here is the result without any pinning pressure

And here is the result for a pinning pressure of 0p01

my residual part looks like this

Real
AdvanceFrictionPressure::computeQpResidual()
{

  std::vector<Real> other_ops(_op_num);
  other_ops = assignOtherOps();

  Real SumEtaj = 0.0;
  Real SumEtaj_G = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumEtaj += (*_vals[i])[_qp] * (*_vals[i])[_qp];
    SumEtaj_G += (*_vals[i])[_qp];
  }

  Real g = _g[_qp];
  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];
  }

  Real total = (_test[_i][_qp] * _mu[_qp] *
                (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
               (_grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp]);

  Real friction = _test[_i][_qp] * 3 * g * _u[_qp] * SumEtaj_G;

  Real k = 10000;
  Real sgn = std::tanh(k * total);
  if (total > 0)
    sgn = 1;
  if (total < 0)
    sgn = -1;
  if (total == 0)
    sgn = 0;
  Real sign = (total > 0) - (total < 0);

  Real abs_total = 100 * std::abs(total);
  Real abs_friction = 100 * std::abs(friction);

  Real a = (abs_total - abs_friction);
  Real smooth_condition = 1.0 / (1.0 + std::exp(-k * (a)));


  Real term1 = smooth_condition * (total - sgn * friction);
  Real L = _L[_qp];

  return L * (term1);
}

IS there something I am missing in the kernel. Because it seems the pinning pressure is much less significant in FEM model (MOOSE) as compared to analytical or even FDM model.
Looking forward to hearing from you.
Thank you
Ashish

[GiudGiud]

hello

is the spatial mesh converged?

[tonkmr]

@ashishdhole, I have been able to get better comparisons between the analytical model for a shrinking circular grain. If decreasing the mesh size stops improving it, you should use a smaller time step size and see if that helps.

I think the problem here is that you are equating the weak form with the strong form. The weak form allows us to solve an equation across the mesh in a form equivalent to the strong form, but locally the weak and strong forms are not identical or interchangeable. The expression used to determine the correct expression for P_f' is based on the strong form of f_i', not the weak form. Basically, by using the expressions with the test function and divergence theorem applied you are changing those values from what they should be.

I recommend that you use ACGrGrInterface and ACInterface and make a new kernel. In that new kernel, you will have to calculate the value for f_i' using the expression from the paper, including the Laplacian term:

You can calculate the Laplacian using _second_u[_qp].tr(), where you initialize _second_u in the constructor using _second_u(second()). See CHInterfaceBase.h for an example. This will mean you won't be able to use first order elements, since their Laplacian is zero.

[ashishdhole]

Thank you Michael @tonkmr

Thank you for this valuable suggestion.
so you are saying that the f_i should look like

"(_mu[_qp] * (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
 (_kappa[_qp] * _second_u[_qp].tr())"

instead of

"(_test[_i][_qp] * _mu[_qp] *
                (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
               (_grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp])"

[tonkmr]

Exactly.

[ashishdhole]

One more question (coding part)
so my residual now looks like

Real
AdvanceFrictionPressure::computeQpResidual()
{

  std::vector<Real> other_ops(_op_num);
  other_ops = assignOtherOps();

  Real SumEtaj = 0.0;
  Real SumEtaj_G = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumEtaj += (*_vals[i])[_qp] * (*_vals[i])[_qp];
    SumEtaj_G += (*_vals[i])[_qp];
  }

  Real g = _g[_qp];
  
  Real f_i = (_mu[_qp] * (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
             (_kappa[_qp] * _second_u[_qp].tr());
  Real f_p = 3 * g * _u[_qp] * SumEtaj_G;

  Real total = (_test[_i][_qp] * _mu[_qp] *
                (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
               (_grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp]);

  Real friction = 3 * g * _u[_qp] * SumEtaj_G;

  Real k = 10000;
  Real sgn = std::tanh(k * f_i);
  if (f_i > 0)
    sgn = 1;
  if (f_i < 0)
    sgn = -1;
  if (f_i == 0)
    sgn = 0;

  Real abs_total = std::abs(f_i);
  Real abs_friction = std::abs(f_p);

  Real a = (abs_total - abs_friction);
  Real smooth_condition = 1.0 / (1.0 + std::exp(-k * (a)));
  if (smooth_condition > 0.5)
    smooth_condition = 1;
  if (smooth_condition < 0.5)
    smooth_condition = 0;

  Real term1 = smooth_condition * (total - sgn * friction);

  return _L[_qp] * (term1);
}

Here I have smooth_condition which basically makes the residual = 0 if f_p >f_i

in jacobian part I have to differentiate wrt _u[_qp]

in that case the differentiation of f_i and f_p wrt _u should look like

  Real df_i_dopi = (_mu[_qp] * (_phi[_j][_qp] * (3.0 * _u[_qp] * _u[_qp] - 1.0 + 3 * SumEtaj))) +
                   (_kappa[_qp] * _second_phi[_j][_qp].tr());
  Real df_p_dopi = (_phi[_j][_qp] * 3 * g * SumEtaj_G);

is that correct?

[ashishdhole]

Hell Michael @tonkmr @GiudGiud

in my modified code,
to make residual = 0 when f_i <= f_p I am introducing a variable called smooth_condition, which is a function of f_i and f_p,

without that condition, the code works fine even when the pinning pressure (g) != 0,

But on introducing the smooth_condition the simulation does not converge. I tried tightning the tolerances and refining the mesh as well. Can you suggest another way to introduce such condition.

I am very close to the actual code, but lack in the knowledge of how moose handles such conditions. I tried to find if there is an existing kernel which uses such condition, but I could not find any.

following is my modified kernel.

#include "AdvanceFrictionPressure.h"

registerMooseObject("pfubcApp", AdvanceFrictionPressure);

InputParameters
AdvanceFrictionPressure::validParams()
{
  InputParameters params = 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 thekernel ");
  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)
  : Kernel(parameters),
    _L(getMaterialProperty<Real>("mob_name")),
    _kappa(getMaterialProperty<Real>("kappa_name")),
    _g(getMaterialProperty<Real>("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")),
    _second_u(second()),
    _second_phi(secondPhi())
{
}

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

  return other_ops;
}

Real
AdvanceFrictionPressure::computeQpResidual()
{

  std::vector<Real> other_ops(_op_num);
  other_ops = assignOtherOps();

  Real SumEtaj = 0.0;
  Real SumEtaj_G = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumEtaj += (*_vals[i])[_qp] * (*_vals[i])[_qp];
    SumEtaj_G += (*_vals[i])[_qp];
  }

  Real g = _g[_qp];

  Real f_i = (_mu[_qp] * (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
             (_kappa[_qp] * _second_u[_qp].tr());
  Real f_p = 3 * g * _u[_qp] * SumEtaj_G;

  Real total = (_test[_i][_qp] * _mu[_qp] *
                (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
               (_grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp]);

  Real friction = _test[_i][_qp] *3 * g * _u[_qp] * SumEtaj_G;

  Real k = 10000;
  Real sgn = std::tanh(k * total);
  if (total > 0)
    sgn = 1;
  if (total < 0)
    sgn = -1;
  if (total == 0)
    sgn = 0;

  Real abs_total = std::abs(f_i);
  Real abs_friction = std::abs(f_p);

  Real a = (abs_total - abs_friction);
  Real smooth_condition = 1.0 / (1.0 + std::exp(-k * (a)));
  if (smooth_condition > 0.5)
    smooth_condition = 1;
  if (smooth_condition < 0.5)
    smooth_condition = 0;

  Real term1 = smooth_condition * (total - sgn * friction);

  return _L[_qp] * (term1);
}

Real
AdvanceFrictionPressure::computeQpJacobian()
{

  std::vector<Real> other_ops(_op_num);
  other_ops = assignOtherOps();

  Real SumEtaj = 0.0;
  Real SumEtaj_G = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumEtaj += (*_vals[i])[_qp] * (*_vals[i])[_qp];
    SumEtaj_G += (*_vals[i])[_qp];
  }

  Real g = _g[_qp];

  Real f_i = (_mu[_qp] * (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
             (_kappa[_qp] * _second_u[_qp].tr());
  Real f_p = 3 * g * _u[_qp] * SumEtaj_G;

  Real df_i_dopi = (_mu[_qp] * (_phi[_j][_qp] * (3.0 * _u[_qp] * _u[_qp] - 1.0 + 3 * SumEtaj))) +
                   (_kappa[_qp] * _second_phi[_j][_qp].tr());
  Real df_p_dopi = (_phi[_j][_qp] * 3 * g * SumEtaj_G);

  Real total = (_test[_i][_qp] * _mu[_qp] *
                (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
               (_grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp]);

  Real friction = _test[_i][_qp] * 3 * g * _u[_qp] * SumEtaj_G;

  Real k = 10000;
  Real sgn = std::tanh(k * total);
  if (total > 0)
    sgn = 1;
  if (total < 0)
    sgn = -1;
  if (total == 0)
    sgn = 0;

  Real abs_total = std::abs(f_i);
  Real abs_friction = std::abs(f_p);

  Real a = (abs_total - abs_friction);
  Real smooth_condition = 1.0 / (1.0 + std::exp(-k * (a)));
  if (smooth_condition > 0.5)
    smooth_condition = 1;
  if (smooth_condition < 0.5)
    smooth_condition = 0;

  Real dtotal_dopi =
      (_test[_i][_qp] *
       (_mu[_qp] * (_phi[_j][_qp] * (3.0 * _u[_qp] * _u[_qp] - 1.0 + 3 * SumEtaj)))) +
      (_grad_phi[_j][_qp] * _kappa[_qp] * _grad_test[_i][_qp]);
  Real dfriction_dopi = ((_test[_i][_qp] *_phi[_j][_qp] * 3 * g * SumEtaj_G));

  Real dsgn_dopi = k * (1.0 - sgn * sgn) * (df_i_dopi);

  Real sign_total = std::tanh(k * df_i_dopi);
  Real sign_friction = std::tanh(k * df_p_dopi);

  Real dsmooth_condition_dopi = k * smooth_condition * (1.0 - smooth_condition) *
                                (sign_total * df_i_dopi - sign_friction * df_p_dopi);

  Real term2 = smooth_condition * (dtotal_dopi - dsgn_dopi * friction - sgn * dfriction_dopi) +
               dsmooth_condition_dopi * (total - sgn * friction);

  return _L[_qp] * term2;
}

Real
AdvanceFrictionPressure::computeQpOffDiagJacobian(unsigned int jvar)
{

  std::vector<Real> other_ops(_op_num);
  other_ops = assignOtherOps();

  Real SumEtaj = 0.0;
  Real SumEtaj_G = 0.0;
  for (unsigned int i = 0; i < _op_num; ++i)
  {
    SumEtaj += (*_vals[i])[_qp] * (*_vals[i])[_qp];
    SumEtaj_G += (*_vals[i])[_qp];
  }

  Real g = _g[_qp];

  Real f_i = (_mu[_qp] * (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
             (_kappa[_qp] * _second_u[_qp].tr());
  Real f_p = 3 * g * _u[_qp] * SumEtaj_G;

  Real total = (_test[_i][_qp] * _mu[_qp] *
                (_u[_qp] * _u[_qp] * _u[_qp] - _u[_qp] + 3 * _u[_qp] * SumEtaj)) +
               (_grad_u[_qp] * _kappa[_qp] * _grad_test[_i][_qp]);

  Real friction = _test[_i][_qp] * 3 * g * _u[_qp] * SumEtaj_G;

  Real k = 10000;
  Real sgn = std::tanh(k * f_i);
  if (f_i > 0)
    sgn = 1;
  if (f_i < 0)
    sgn = -1;
  if (f_i == 0)
    sgn = 0;

  Real abs_total = std::abs(f_i);
  Real abs_friction = std::abs(f_p);

  Real a = (abs_total - abs_friction);
  Real smooth_condition = 1.0 / (1.0 + std::exp(-k * (a)));
  if (smooth_condition > 0.5)
    smooth_condition = 1;
  if (smooth_condition < 0.5)
    smooth_condition = 0;

  for (unsigned int i = 0; i < _op_num; ++i)
    if (jvar == _vals_var[i])
    {
      Real dSumOPj = 2.0 * (*_vals[i])[_qp] * _phi[_j][_qp];
      Real dSumOPGj = _phi[_j][_qp];
      Real dtotal_dopj = (_test[_i][_qp] * ((_mu[_qp] * 3 * _u[_qp] * dSumOPj)));
      Real dfriction_dopj = ((_test[_i][_qp] * 3 * g * _u[_qp] * dSumOPGj));

      Real df_i_dopj = (_mu[_qp] * 3 * _u[_qp] * dSumOPj);
      Real df_p_dopj = (3 * g * _u[_qp] * dSumOPGj);

      Real dsgn_dopj = k * (1.0 - sgn * sgn) * (df_i_dopj);

      Real sign_total = std::tanh(k * df_i_dopj);
      Real sign_friction = std::tanh(k * df_p_dopj);

      Real dsmooth_condition_dopj = k * smooth_condition * (1.0 - smooth_condition) * 
                                    (sign_total * df_i_dopj - sign_friction * df_p_dopj);

      Real term3 = smooth_condition * (dtotal_dopj - dsgn_dopj * friction - sgn * dfriction_dopj) +
                   dsmooth_condition_dopj * (total - sgn * friction);

      return _L[_qp] * term3;
    }

  return 0;
}

you can see I have used the step function for sgn and smooth_condition, which could cause non convergence, but, without them also the problem does not converge when g != 0. with g = 0, the code works fine.

even if it converges, the time step becomes very small (1e-7). I tried to reduce the tolerance to 1e-9. The pinning pressure I have used is so small that the grain growth should happen.

Any suggestion will be highly appreciated.
Thank you