RheologyCalculator.jl: automatically solving arbitrary constitutive relationships

Much of the challenges in simulating geoscientific processes stems from the (often) non-linear nature and complexity of the constitutive equations. Rheological behavior is often modeled as combinations of individual elements (such as linear elastic springs, viscous dashpots, and plastic sliders), each representing a simple (nonlinear) constitutive relation. These elements are assembled in series and/or parallel to form more complex rheological systems. Solving the resulting constitutive equations typically requires iterative methods, such as Newton-type schemes, which rely on the assembly of local Jacobian matrices.

In most existing codes, constitutive relationships and their corresponding Jacobians are hard-coded, making the implementation of new rheologies a challenging task, as it requires deep modifications to the codebase. Here, we present RheologyCalculator.jl, a Julia package designed to solve constitutive relationships constructed from arbitrary combinations of rheological elements. The package statically determines the computational graph of the constitutive model at compile time, enabling efficient assembly of local Jacobian matrices and the use of local direct solvers without heap allocations, which is crucial for maintaining high performance when employed in large-scale simulations, as the system of equations must be solved in a large number of grid points. RheologyCalculator.jl is also fully differentiable, allowing it to, for example, compute consistent tangent operators or to be differentiated through entire forward simulations for gradient-based analyses. The package is designed for seamless integration into existing geoscientific modeling frameworks, providing a flexible and efficient approach for implementing and experimenting with complex rheological models without extensive code modifications.