How can I add external C++ libraries to use in MOOSE?

[richmondodufisan]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

First I want to know generally how I would add an external library and add it to the Makefile so I can compile. Whether autodiff (what I'm asking about) or any other library in general. Also curious if there is a way to write a python script which is called, does the calculations I want, and returns the results to the C++ source file.

Second, I want to know how to add autodiff (https://github.com/autodiff/autodiff) and wondering if it necessary since MOOSE already has some kind of AD system that is supposed to do the same thing? I'm not sure. In my custom Material file I am currently calculating derivatives manually with finite differencing, for example:

// Compute the numerical derivative of a scalar function w.r.t a tensor_ij

RankTwoTensor ComputeStrainEnergyNeoHookeanNearlyIncompressible::compute_dWdC(const Real &mu_0, const RankTwoTensor &C) 
{
	// Initialize perturbation
	Real epsilon = 1e-6;
	
	// Initialize derivative
	RankTwoTensor dWdC;
	
    // Compute the scalar function for the original tensor
    Real W = computeStrainEnergy(mu_0, C);

    // Loop over all components of the tensor to compute the derivative
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            // Create a copy of the tensor to perturb
            RankTwoTensor C_perturbed = C;

            // Perturb the (i, j) component by epsilon
            C_perturbed(i, j) += epsilon;

            // Evaluate the scalar function at the perturbed tensor
            Real W_perturbed = computeStrainEnergy(mu_0, C_perturbed);

            // Compute the finite difference derivative and store it in the dfdT tensor
            dWdC(i, j) = (W_perturbed - W) / epsilon;
        }
    }

    // Return the tensor of derivatives
    return dWdC;
}

And I am getting good results, but I am wondering if the convergence and accuracy is affected by this because it solves a bit slowly and my choice of epsilon is problem-dependent. I want to instead add autodiff (https://github.com/autodiff/autodiff) do something like (I used chatGPT for this sample but I want to do something similar):

RankTwoTensor ComputeStrainEnergyNeoHookeanNearlyIncompressible::compute_dWdC(const Real &mu_0, const RankTwoTensor &C)
{
    using namespace autodiff;
    using dual = dual2nd;

    Eigen::Matrix<dual, 3, 3> C_dual;
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            C_dual(i, j) = C(i, j);

    auto W = [&](const Eigen::Matrix<dual, 3, 3> & C_dual) {
        dual J = std::pow(C_dual.determinant(), 1.0 / 2.0);
        Eigen::Matrix<dual, 3, 3> C_iso = std::pow(J, -2.0 / 3.0) * C_dual;
        dual I_1 = C_iso.trace();
        return (mu_0 / 2.0) * (I_1 - 3);
    };

    Eigen::Matrix<dual, 3, 3> dWdC_dual = gradient(W, wrt(C_dual), at(C_dual));
    RankTwoTensor dWdC;
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            dWdC(i, j) = val(dWdC_dual(i, j));  // Extract the value from dual

    return dWdC;
}

How can I do this? I have already tried https://mooseframework.inl.gov/source/materials/DerivativeMaterialInterface.html and it didn't work to calculate derivatives of scalars wrt tensors

[GiudGiud]

Hello

First I want to know generally how I would add an external library and add it to the Makefile so I can compile.

If the library is already compiled on your machine, you'll want to add to several environment variables in the Makefile.
However, I dont think you need to do that here.

We have automatic differentiation in moose using MetaPhyscL through libmesh. Every single object in MOOSE prefaced with AD already uses this. And all you need to do to use it is:

• inherit from the ADMaterial class for your custom material
• everywhere you have Real, use ADReal, same for RealVectorValue, RealTensorValue etc

[richmondodufisan]

it compiled but didn't converge unfortunately. I added a cout to see the size of W.derivatives() and it was extremely large, 12816. I don't think the way I set it up is correct.

[GiudGiud]

ADReal W = computeStrainEnergy(mu_0, C_seeded);

mu_0 does not contain any derivatives with regards to nonlinear variables right?
and neither does any term that comes in computeStrainEnergy right?

We cant get the regular use of AD and our local use mixed up

[richmondodufisan]

Nope, mu_0 is just a constant ADReal. And the function is here, could be because of all the matrix operations?

ADReal ComputeStrainEnergyNeoHookeanNearlyIncompressible::computeStrainEnergy(const ADReal &mu_0, const ADRankTwoTensor &C)
{
	ADReal J = std::pow(C.det(), 1.0 / 2.0);
	
	// Isochoric Right Cauchy-Green deformation tensor
	ADRankTwoTensor C_iso = std::pow(J, -2.0 / 3.0) * C;
	
	// First invariant, equivalent to lambda_1^2 + lambda_2^2 + lambda_3^2
	ADReal I_1 = C_iso.trace();
	
	ADReal strain_energy = (mu_0/2.0) * (I_1 - 3);
	
	return strain_energy;
}

[GiudGiud]

this should be fine. However the "det()" operation likely involves every single component, making every component depend on each other

[richmondodufisan]

I see. I guess numerical differentiation is the best option here then, I just have to toggle the perturbation value each time. thanks for your help!