Calculation of pk2_jacobian in moose

[manikumd22]

Hi @GiudGiud ,

In MOOSE, when calculating the PK2 Jacobian (C=∂S/∂E), does the framework use automatic differentiation? Specifically, if I directly compute the Second Piola-Kirchhoff stress (S) in my model, will MOOSE automatically differentiate S with respect to the Green-Lagrange strain (E) to compute the Jacobian (C)?

[maxnezdyur]

No that system does not use AD. It does have perfect handcoded jacobians. I am not sure what you mean by "directly computing" S. If you changed how you compute S then I think you will need to fix the pk2 jacobian too.

[reverendbedford]

The material objects will calculate the strain, all you have to do is calculate the stress and derivative.

For PK1 we don't have a "native" material model using this interface directly. But the PK2 actually derives from it, so you can see what you need to set here. You "get" the deformation gradient and need to provide the PK1 stress and the derivative of PK1 with respect to the deformation gradient.

For PK2 there is a complete example here. You "get" the Green-Lagrange strain and have to provide the PK2 stress and the derivative of PK2 with respect to the deformation gradient.

[manikumd22]

Hi @GiudGiud

As for the Pk1 kernel
this is formulation

If I calculate _pk1_stress[_qp], then where it is taking derivative?

if (_large_kinematics)
{
_cauchy_stress[_qp] = _pk1_stress[_qp] * _F[_qp].transpose() / _F[_qp].det();

auto f = _inv_df[_qp].inverse();
usingTensorIndices(i_, j_, k_, l_);
_cauchy_jacobian[_qp] = _cauchy_stress[_qp].times<i_, l_, j_, k_>(f) -
                        _cauchy_stress[_qp].times<i_, j_, k_, l_>(f.transpose());
_cauchy_jacobian[_qp] +=
    _pk1_jacobian[_qp].tripleProductJkl(_F[_qp], f.transpose(), _F[_qp]) / _F[_qp].det();

}

[reverendbedford]

You don't worry about the kernel, you just need to worry about the material. The rest of the system takes care of calculating the derivative the kernel needs.

So if you want to inherit from the ComputeLagrangianStressPK1 object all your material needs to do is set:

• _pk1_stress[_qp] to the PK1 stress
• _pk1_jacobian[_qp] to the derivative of the PK1 stress with respect to the deformation gradient.

This will happen in your material object, which again should derive from ComputeLagrangianStressPK1

[manikumd22]

Hi @reverendbedford

I have calculated the _pk1stress and perturbed _pk1 stress

this is how I will calculate the _pk1jacobian[_qp], right?

            RankTwoTensor original_F = _F[_qp]; // Original deformation gradient
            RankTwoTensor original_P = _pk1_stress[_qp]; // Original stress tensor
            RankTwoTensor perturbed_P = _perturbed_P[_qp]; // Original stress tensor
            RankFourTensor C; // Fourth-order elasticity tensor to store results

            // Loop over all components of F
            for (unsigned int k = 0; k < 2; ++k) // Assume 2x2 for 2D (adjust for 3D as needed)
            {
                for (unsigned int l = 0; l < 2; ++l)
                {
                    // Perturb the (k, l) component of F
                    RankTwoTensor perturbed_F = original_F;
                    perturbed_F(k, l) += delta;

                    // Compute finite difference derivative
                    for (unsigned int i = 0; i < 2; ++i) // Rows of P
                    {
                        for (unsigned int j = 0; j < 2; ++j) // Columns of P
                        {
                            // C_{ijkl} = (P_{ij}(perturbed) - P_{ij}(original)) / delta
                            C(i, j, k, l) = (perturbed_P(i, j) - original_P(i, j)) / (delta*10e9);
                        }
                    }
                }
            }

            _pk1_jacobian[_qp] = C;

[reverendbedford]

Uhh, this is something you'll probably have to debug in detail yourself, but a couple of things stand out:

1. All MOOSE (solid mechanics) materials should be setup in 3D, even if your simulation is 2D or 1D.
2. The idea of perturbing the deformation gradient is right, but where are you actually calling your model again to get the perturbed PK1 stress?
3. Not sure where the magic 10e9 is coming from here.