Finding computeQpOffDiagJacobian for sink and recombination #27660

jackeggemeyer · 2024-05-17T20:41:18Z

jackeggemeyer
May 17, 2024

I am trying to add sink and recombination to the Multi-phase poly-crystalline Allen-Cahn Kernel. The top equation is the equation already used by MOOSE. The second equation is from: https://www.sciencedirect.com/science/article/pii/S0927025622004839#! and it adds vacancy production, vacancy-interstial recombination, and sink terms to the grand potential equation in terms of chemical potential.

From SuspectabilityTimeDirivative, matdiffusion, maskedbodyforce(deals with h_m*s_v^o term), and coupledswitchingtimedirivative I was able to recreate the equation from image 1:

(_test[_i][_qp] * _u_dot[_qp]) * _Chi[_qp]; = _D[_qp] * _grad_v[_qp] * _grad_test[_i][_qp]; +(-_test[_i][_qp] * _scale * _postprocessor * _function.value(_t, _q_point[_qp]);) * _mask[_qp]; - ((*_prop_dhjdetai[n])[_qp] * (*_prop_Fj[n])[_qp];)) * (_test[_i][_qp] * _v_dot[_qp];)

I am creating a new kernal that adds the missing term:

h_m * K_iv * rho * rho1 - h_m * k_vs * rho     (the rho1+rho2 is just rho1)

Here is the residual in moose i created:

 ((*_prop_hm[n])[_qp] * (*_prop_kiv[n]) (*_prop_Rho[n])[_qp] * (*_prop_Rhos[n])[_qp] + (*_prop_hm[n])[_qp] * (*_prop_kxs[n]) * (*_prop_Rho[n])[_qp];)

I copied below my .c file, I am struggling to find documentation on finding my off diagonal Jacobian equation. I have been following the steps provided by the INL: https://mooseframework.inl.gov/getting_started/examples_and_tutorials/tutorial01_app_development/step04_weak_form.html#procedure
I am still struggling to confirm what my equations should look like. Any tips or documentation on how to create your offdiag would be awesome. I attached the .c and .h in text files and pasted the .c below! Also this is my first kernel ever so if you see any other issues please let me know :).

p.s. Am I missing an override anywhere?

#include "GrandPotentialSink.h"

registerMooseObject("PhaseFieldApp", GrandPotentialSink);

InputParameters
GrandPotentialSink::validParams()
{
  InputParameters params = CoupledTimeDerivative::validParams();
  params.addClassDescription("Adds sink and recombination to grand potential multi-phase model.");
  params.addRequiredParam<MaterialPropertyName>("Rho", "Material property for Rho.");
  params.addRequiredParam<MaterialPropertyName>("Rhos", "Material property for Rhos.");
  params.addRequiredParam<MaterialPropertyName>("hm", "Material property for the switching function h.");
  params.addRequiredParam<MaterialPropertyName>("kiv", "Material property for the rate of vacancy-interstitial recombination.");
  params.addRequiredParam<MaterialPropertyName>("kxs", "Material property for the product of the rate of vacancy absorption at sinks and sink density.");
  params.addCoupledVar("v", "Order parameter variable.");
  return params;
}

GrandPotentialSink::GrandPotentialSink(const InputParameters & parameters)
  : DerivativeMaterialInterface<JvarMapKernelInterface<CoupledTimeDerivative>>(parameters),
    _v_name(coupledName("v")),
    _prop_Rho(getMaterialProperty<Real>("Rho")),
    _prop_dRhodv(getMaterialPropertyDerivative<Real>("Rho", _v_name)),
    _prop_Rhos(getMaterialProperty<Real>("Rhos")),
    _prop_dRhosdv(getMaterialPropertyDerivative<Real>("Rhos", _v_name)),
    _prop_hm(getMaterialProperty<Real>("hm")),
    _prop_dhmdetai(getMaterialPropertyDerivative<Real>("hm", _v_name)),
    _prop_d2hmdetai2(getMaterialPropertyDerivative<Real>("hm", _v_name, _v_name)),
    _prop_kiv(getMaterialProperty<Real>("kiv")),
    _prop_dkivdetai(getMaterialPropertyDerivative<Real>("kiv", _v_name)),
    _prop_d2kivdetai2(getMaterialPropertyDerivative<Real>("kiv", _v_name, _v_name)),
    _prop_kxs(getMaterialProperty<Real>("kxs")),
    _prop_dkxsdetai(getMaterialPropertyDerivative<Real>("kxs", _v_name)),
    _prop_d2kxsdetai2(getMaterialPropertyDerivative<Real>("kxs", _v_name, _v_name))
{
}

void
GrandPotentialSink::initialSetup()
{
  validateNonlinearCoupling<Real>("Rho");
  validateNonlinearCoupling<Real>("hm");
  validateNonlinearCoupling<Real>("Rhos");
  validateNonlinearCoupling<Real>("kiv");
  validateNonlinearCoupling<Real>("kxs");
}

Real
GrandPotentialSink::computeQpResidual()
{
  return _prop_hm[_qp] * _prop_kiv[_qp] * _prop_Rho[_qp] * _prop_Rhos[_qp] + _prop_hm[_qp] * _prop_kxs[_qp] * _prop_Rho[_qp];
}

Real
GrandPotentialSink::computeQpJacobian()
{
  return _prop_hm[_qp] * _prop_kiv[_qp] *_prop_dRhodv[_qp] * _prop_Rhos[_qp] * _phi[_j][_qp] + 
  _prop_hm[_qp] * _prop_kiv[_qp] * _prop_Rho[_qp] * _prop_dRhosdv[_qp] * _phi[_j][_qp] + 
  _prop_hm[_qp] * _prop_kxs[_qp] * _prop_dRhodv[_qp] * _phi[_j][_qp]
}

Real
GrandPotentialSink::computeQpOffDiagJacobian(unsigned int jvar)
{
  return _prop_dhmdetai[_qp] * _prop_kiv[_qp] * _prop_Rho[_qp] * _prop_Rhos[_qp] + _prop_dhmdetai[_qp] * _prop_kxs[_qp] * _prop_Rho[_qp];
  
}

PotentialSink.txt

PotentialSinkhfile.txt

Answered by GiudGiud

May 29, 2024

Hello

It s missing this header : JvarMapKernelInterface.h

Does it exist in your version of moose?
Can you find it in either the framework or the modules headers?

Guillaume

View full answer

GiudGiud · 2024-05-19T23:27:07Z

GiudGiud
May 19, 2024
Collaborator

Hello

((*_prop_hm[n])[_qp] * (*_prop_kiv[n]) (*_prop_Rho[n])[_qp] * (*_prop_Rhos[n])[_qp] + (*_prop_hm[n])[_qp] * (*_prop_kxs[n]) * (*_prop_Rho[n])[_qp];)

what is n here? I see _qp (which seems correct) in your kernel

For the Jacobian:

do only rho and rho_s depend on variable v? h_m, k_iv and k_vs do not?

For the off-diagonal

does only h-m depend on eta? rho,rho_s and k_vs do not?
I would expect an if-condition there, checking for jvar = index for eta

p.s. Am I missing an override anywhere?

All of these should be override. The nice thing about override is that if you made a mistake in the signature (for example the return type), it wont compile thanks to override

  virtual void initialSetup();

protected:
  virtual Real computeQpResidual();
  virtual Real computeQpJacobian();
  virtual Real computeQpOffDiagJacobian(unsigned int jvar);

NB: do not forget to check your Jacobian with the Jacobian testing tools in MOOSE. See JacobianTest on the website, or troubleshooting failed solves

5 replies

jackeggemeyer May 20, 2024
Author

Thank you so much for your help, the [n] was a mistake as I used CoupledSwitchingTimeDirivative as a template and it inputted multiple values for each and I didnt notice I was pasting an old version of the equation. Here is the updated one:

Real
GrandPotentialSink::computeQpResidual()
{
  return _prop_hm[_qp] * _prop_kiv[_qp] * _prop_Rho[_qp] * _prop_Rhos[_qp] + _prop_hm[_qp] * _prop_kxs[_qp] * _prop_Rho[_qp];
}

For the Jacobian:
Yes, only rho and rho_s depend on variable v and they are DerivativeParsedMaterials, here is what they look like in my input file:

[./rhovbub]
    type = DerivativeParsedMaterial
    args = 'wv'
    f_name = rhovbub
    material_property_names = 'Va kvbub cvbubeq'
    function = 'wv/Va^2/kvbub + cvbubeq/Va'
    derivative_order = 2
  [../]
  [./rhovmatrix]
    type = DerivativeParsedMaterial
    args = 'wv'
    f_name = rhovmatrix
    material_property_names = 'Va kvmatrix cvmatrixeq'
    function = 'wv/Va^2/kvmatrix + cvmatrixeq/Va'
    derivative_order = 2
  [../]
  [./rhov]
    type = DerivativeParsedMaterial
    args = 'wv'
    f_name = rhov
    material_property_names = 'hm hb rhovmatrix(wv) rhovbub(wv)'
    function = 'hm * rhovmatrix + hb * rhovbub'
    derivative_order = 2
    output_properties = 'rhov'
    outputs = exodus
  [../]

K_iv and K_vs are thankfully constants and not dependent on either eta or v.

h_m is a SwitchingFunctionMultiPhaseMaterial and is dependant on eta.

  [./hm]
    type = SwitchingFunctionMultiPhaseMaterial
    h_name = hm
    all_etas = 'etab0 etam0'
    phase_etas = 'etam0'
  [../]

If those are the only variables that effect it, the above Jacobian should be correct right?

Yes, eta is the only variable effecting h_m as seen above. I added the if statement but I am still confused on what the offdiagjaacobina should look like. I have been following both beginner moose tutorials and the step were they describe going from the weak form too the offdiagjacobian still confuses me. Here is what I have so far:

GrandPotentialSink::computeQpOffDiagJacobian(unsigned int jvar)
{
  const unsigned int cvar = mapJvarToCvar(jvar);

  if (jvar == _v_var)
  {
    Real sum = 0.0;

      sum = _prop_dhmdeta[_qp] * _prop_kiv[_qp] * _prop_Rho[_qp] * _prop_Rhos[_qp] + _prop_dhmdetai[_qp] * _prop_kxs[_qp] * _prop_Rho[_qp];

    return _phi[_j][_qp] * sum * _test[_i][_qp];
  }

  Real sum = 0.0;
    sum = _prop_dhmdarg[_qp] * _prop_kiv[_qp] * _prop_Rho[_qp] * _prop_Rhos[_qp] + _prop_dhmdargi[_qp] * _prop_kxs[_qp] * _prop_Rho[_qp];

  return sum * _phi[_j][_qp];
  
}

Thank you so much for the help, let me know if you need anything else.

And thank you so much for pointing out the Jacobian Testing Tool thats extremely helpful

GiudGiud May 21, 2024
Collaborator

what are _prop_dhmdarg and _prop_dhmdargi?
which variables are they the derivative with regards to?
and why is jvar not being used to limit this 'return this sum' to the indexes for that variable?

also how do you set _v_var ?

jackeggemeyer May 28, 2024
Author

Apologies, i forgot to delete some of the previous variable names. Here is my updated file:

#include "GrandPotentialSink.h"

registerMooseObject("PhaseFieldApp", GrandPotentialSink);

InputParameters
GrandPotentialSink::validParams()
{
  InputParameters params = CoupledTimeDerivative::validParams();
  params.addClassDescription("Adds sink and recombination to grand potential multi-phase model.");
  params.addRequiredParam<MaterialPropertyName>("Rho", "Material property for Rho.");
  params.addRequiredParam<MaterialPropertyName>("Rhos", "Material property for Rhos.");
  params.addRequiredParam<MaterialPropertyName>("hm", "Material property for the switching function h.");
  params.addRequiredParam<MaterialPropertyName>("kiv", "Material property for the rate of vacancy-interstitial recombination.");
  params.addRequiredParam<MaterialPropertyName>("kxs", "Material property for the product of the rate of vacancy absorption at sinks and sink density.");
  params.addCoupledVar("v", "Order parameter variable.");
  return params;
}

GrandPotentialSink::GrandPotentialSink(const InputParameters & parameters)
  : DerivativeMaterialInterface<JvarMapKernelInterface<CoupledTimeDerivative>>(parameters),
    _v_var(coupledValue("v")),
    _v_var_num(coupled("v")),
    Rho(getMaterialProperty<Real>("Rho")),
    dRhodv(getMaterialPropertyDerivative<Real>("Rho", coupled("v"))),
    Rhos(getMaterialProperty<Real>("Rhos")),
    dRhosdv(getMaterialPropertyDerivative<Real>("Rhos", coupled("v"))),
    hm(getMaterialProperty<Real>("hm")),
    dhmdeta(getMaterialPropertyDerivative<Real>("hm", coupled("v"))),
    d2hmdeta2(getMaterialPropertyDerivative<Real>("hm", coupled("v"), coupled("v"))),
    kiv(getMaterialProperty<Real>("kiv")),
    kxs(getMaterialProperty<Real>("kxs"))
{
}

void
GrandPotentialSink::initialSetup()
{
  validateNonlinearCoupling<Real>("Rho");
  validateNonlinearCoupling<Real>("hm");
  validateNonlinearCoupling<Real>("Rhos");
  validateNonlinearCoupling<Real>("kiv");
  validateNonlinearCoupling<Real>("kxs");
}

Real
GrandPotentialSink::computeQpResidual()
{
  return hm[_qp] * kiv[_qp] * Rho[_qp] * Rhos[_qp] + hm[_qp] * kxs[_qp] * Rho[_qp];
}

Real
GrandPotentialSink::computeQpJacobian()
{
  return hm[_qp] * kiv[_qp] * dRhodv[_qp] * Rhos[_qp] * _phi[_j][_qp] + 
         hm[_qp] * kiv[_qp] * Rho[_qp] * dRhosdv[_qp] * _phi[_j][_qp] + 
         hm[_qp] * kxs[_qp] * dRhodv[_qp] * _phi[_j][_qp];
}

Real
GrandPotentialSink::computeQpOffDiagJacobian(unsigned int jvar)
{
  if (jvar == _v_var_num)
  {
    Real sum = dhmdeta[_qp] * kiv[_qp] * Rho[_qp] * Rhos[_qp] +
               hm[_qp] * kiv[_qp] * dRhodv[_qp] * Rhos[_qp] +
               hm[_qp] * kiv[_qp] * Rho[_qp] * dRhosdv[_qp] +
               dhmdeta[_qp] * kxs[_qp] * Rho[_qp] +
               hm[_qp] * kxs[_qp] * dRhodv[_qp];

    return _phi[_j][_qp] * sum * _test[_i][_qp];
  }

  return 0.0;
}

#pragma once

#include "CoupledTimeDerivative.h"
#include "JvarMapKernelInterface.h"
#include "DerivativeMaterialInterface.h"

class GrandPotentialSink : public DerivativeMaterialInterface<JvarMapKernelInterface<CoupledTimeDerivative>>
{
public:
  static InputParameters validParams();
  GrandPotentialSink(const InputParameters & parameters);

protected:
  virtual void initialSetup() override;
  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobian() override;
  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;

private:
  const VariableValue & _v_var;
  const unsigned int _v_var_num;
  const MaterialProperty<Real> & Rho;
  const MaterialPropertyDerivative<Real> & dRhodv;
  const MaterialProperty<Real> & Rhos;
  const MaterialPropertyDerivative<Real> & dRhosdv;
  const MaterialProperty<Real> & hm;
  const MaterialPropertyDerivative<Real> & dhmdeta;
  const MaterialPropertyDerivative<Real> & d2hmdeta2;
  const MaterialProperty<Real> & kiv;
  const MaterialProperty<Real> & kxs;
};

I don't understand the jvar if statements that well, I understand they have to do with your coupled variables so if I am misusing it apologies. I properly set the _v_var variable and I am now getting this error:

(moose) jeggeme@jeggeme-OptiPlex-Tower-7010:~/projects/tarpon$ make -j 6
Checking if header needs updating: /home/jeggeme/projects/moose/framework/include/base/MooseRevision.h...
Compiling C++ (in opt mode) /home/jeggeme/projects/tarpon/build/unity_src/kernels_Unity.C...
In file included from /home/jeggeme/projects/tarpon/src/kernels/GrandPotentialSink.C:1,
                 from /home/jeggeme/projects/tarpon/build/unity_src/kernels_Unity.C:2:
/home/jeggeme/projects/tarpon/build/header_symlinks/GrandPotentialSink.h:4:10: fatal error: JvarMapKernelInterface.h: No such file or directory
    4 | #include "JvarMapKernelInterface.h"
      |          ^~~~~~~~~~~~~~~~~~~~~~~~~~
compilation terminated.
make: *** [/home/jeggeme/projects/moose/framework/build.mk:151: /home/jeggeme/projects/tarpon/build/unity_src/kernels_Unity.x86_64-conda-linux-gnu.opt.lo] Error 1
(moose) jeggeme@jeggeme-OptiPlex-Tower-7010:~/projects/tarpon$

GiudGiud May 29, 2024
Collaborator

Hello

It s missing this header : JvarMapKernelInterface.h

Does it exist in your version of moose?
Can you find it in either the framework or the modules headers?

Guillaume

Answer selected by jackeggemeyer

jackeggemeyer May 30, 2024
Author

Yup that was it, Should have checked that. Thank you so much.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Finding computeQpOffDiagJacobian for sink and recombination #27660

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 5 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Select a reply

Finding computeQpOffDiagJacobian for sink and recombination #27660

jackeggemeyer May 17, 2024

Replies: 1 comment · 5 replies

GiudGiud May 19, 2024 Collaborator

jackeggemeyer May 20, 2024 Author

GiudGiud May 21, 2024 Collaborator

jackeggemeyer May 28, 2024 Author

GiudGiud May 29, 2024 Collaborator

jackeggemeyer May 30, 2024 Author

jackeggemeyer
May 17, 2024

Replies: 1 comment 5 replies

GiudGiud
May 19, 2024
Collaborator

jackeggemeyer May 20, 2024
Author

GiudGiud May 21, 2024
Collaborator

jackeggemeyer May 28, 2024
Author

GiudGiud May 29, 2024
Collaborator

jackeggemeyer May 30, 2024
Author