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

ashishdhole · 2024-12-04T22:17:34Z

ashishdhole
Dec 4, 2024

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 · 2024-12-04T22:20:39Z

GiudGiud
Dec 4, 2024
Collaborator

hello

is the spatial mesh converged?

25 replies

ashishdhole Dec 8, 2024
Author

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

ashishdhole Dec 13, 2024
Author

Hello @laagesen ,

Can you have a look at this code for me please.
Thank you

tonkmr Dec 16, 2024

@ashishdhole Ifs often lead to convergence problems. Smoother functions that provide similar functionality may provide better performance. The bad convergence may also be due to mistakes in your Jacobian. I recommend trying things like using a small problem and switching to the Finite Difference preconditioner (https://mooseframework.inl.gov/source/preconditioners/FiniteDifferencePreconditioner.html) to see if that improves things. If it does, that tells you that there is a problem with your jacobian. Note you should never us FDP for an actual solution, just for debugging. If there is a problem with your Jacobian, you could use the Jacobian debugger to try to find it (https://mooseframework.inl.gov/help/development/analyze_jacobian.html).

ashishdhole Dec 16, 2024
Author

Hello @tonkmr I tred with finite difference preconditioners but results are almost similar. and the jacobian debugger is also detecting no errors.

I think the problem is the interpolation function

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

The value calculated by moose is way too less than the driving force of the system. because of which the grain boundaries see no resistance at all. And as the simulations is expected to do, where the grain boundary should stop moving once this resistance is more than the driving force, never stops. I could see the difference more than 100 times.

On the other hand in FDM model, this function seems to be having reasonable value.

ashishdhole Dec 17, 2024
Author

I tried the AD version as well so I don't think the jacobian is the issue.

If you can help me with the variable sgn and smooth_function in my kernel.

There are many functions and variable predefined in moose and many of them are not documented. It is hard to use the exact variables that is accepted by moose framework.

like in the above equation 1st and 3rd term is basically ACGrGrPoly and ACInterface. witout 2nd term my code works just like a grain growth kernel provided by moose. But when I am introducing 2nd term, the kernel does not behaves as it should be.

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

value of this term is very small when I am calculating in moose vs when I am calculating in FDM.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 25 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

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

ashishdhole Dec 4, 2024

Replies: 1 comment · 25 replies

GiudGiud Dec 4, 2024 Collaborator

ashishdhole Dec 8, 2024 Author

ashishdhole Dec 13, 2024 Author

tonkmr Dec 16, 2024

ashishdhole Dec 16, 2024 Author

ashishdhole Dec 17, 2024 Author

ashishdhole
Dec 4, 2024

Replies: 1 comment 25 replies

GiudGiud
Dec 4, 2024
Collaborator

ashishdhole Dec 8, 2024
Author

ashishdhole Dec 13, 2024
Author

ashishdhole Dec 16, 2024
Author

ashishdhole Dec 17, 2024
Author