Interpolation of an expression is slow, starting at a certain time step

General
1

Hello everyone,

I’ll give a brief background about myself, if my post is against the rules, then I guess it’ll get removed.

I’m an undergrad student of materials science. I’ve been tasked with characterizing polymer materials (choosing the right models and fitting the parameters), which will be used in a more involved simulation later on.

I’m very new to the field of finite element analysis. I don’t have access to commercial software, so I had to look for open world solutions and as I have some familiarity with Python already, I chose Fenicsx. I’m slowly getting hang of how things work, but still a lot of things are obscure to me.

My material to investigate is an elastomer, so my first choice for a suitable material model is the Bergstrom-boyce model as presented here:

I believe my implementation of the model itself is correct, at least it seems to produce reasonable animations of tensile, compression, buckling etc. tests in Paraview. However, when I try to interpolate the expression for the stress tensor, at a certain timestep during the simulation the performance is crippled, the iteration being two orders of magnitude slower than the previous one, eventually running into JIT timeout errors. I am guessing I am running into some memory or cache problems, but I have no idea where to even begin to investigate the issue.

Again, as my knowledge of Fenicsx is minimal, most of my solutions are total hackjobs, mostly copy pasted from tutorial. The interpolation approach was copied from this tutorial:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/linearelasticity_code.html

I am posting my code in its entirety, the first half of the script defines the material model and works as intended (I hope), the issue comes specifically from the line:

interpolate_expression(σ)

towards the end of the script, commenting it out resolves the issue of performance woes, but as the stress values are what I’m after, I can’t skip this step. :sweat_smile:

Again, as I’m very new to Fenics and this is my first implementation, I’d be happy to receive feedback if my approach towards any other part of the script is bad/inefficient/faulty.

# imports
from dolfinx import mesh
import numpy as np
from mpi4py import MPI
from time import time
from ufl import (
    Measure,
    TestFunction,
    TrialFunction,
    as_tensor,
    indices,
    Identity,
    grad,
    det,
    tr,
    sqrt,
    dev,
    inv,
    derivative,
)

from dolfinx.fem import (
    VectorFunctionSpace,
    TensorFunctionSpace,
    Function,
    Expression,
    Constant,
    dirichletbc,
    locate_dofs_topological,
)
from dataclasses import dataclass
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx.io import XDMFFile

# Material_BB class definition
@dataclass
class Material_BB:
    """Class for defining parameters for a material employing the Bergstrom-Boyce model"""

    μ: float
    λ_lock: float
    κ: float  # bulk modulus
    s: float  # ratio μB/μA
    ξ: float  # perturbation
    C: float  #
    τ_base: float
    m: float
    τ_cut: float  # cut-off stress below which no flow will occur
    ρ: float  # mass density


def inv_lan(x: float) -> float:
    """Calculates the Padé approximation of the inverse Langevin function.

    Args:
        x (float): scalar argument
    """
    return x * (3 - x**2) / (1 - x**2)


def ramp(x: float) -> float:
    """Calculates the ramp function R.

    Args:
        x (float): scalar argument
    """
    return (x + np.abs(x)) / 2


def get_sigma_A(F, mat: Material_BB):
    J = det(F)
    b = F * F.T
    b_star = J ** (-2 / 3) * b
    λ_starbar = sqrt(tr(b_star) / 3.0)

    σ_A = (mat.μ * inv_lan(λ_starbar / mat.λ_lock) * dev(b_star)) / (
        J * λ_starbar * inv_lan(1 / mat.λ_lock)
    ) + mat.κ * (J - 1) * δ
    return σ_A


def get_sigma_B(F_e: Expression, F_v: Expression, mat: Material_BB):
    J_e = det(F_e)
    b_e = F_e * F_e.T
    b_e_star = J_e ** (-2 / 3) * b_e
    λ_e_starbar = sqrt(tr(b_e_star) / 3.0)

    σ_B = (mat.s * mat.μ * inv_lan(λ_e_starbar / mat.λ_lock) * dev(b_e_star)) / (
        J_e * λ_e_starbar * inv_lan(1 / mat.λ_lock)
    ) + mat.κ * (J_e - 1) * δ

    return σ_B


def get_delta_F_v(F_v, F_e, σ_B, mat: Material_BB, dt):
    J_v = det(F_v)
    b_v = F_v * F_v.T
    b_v_star = J_v ** (-2 / 3) * b_v
    λ_v_starbar = sqrt(tr(b_v_star) / 3.0)
    dev_σ_B = dev(σ_B)
    τ = tr(dev_σ_B.T * dev_σ_B)  # Frobenius norm

    γ_dot = (λ_v_starbar - 1 + mat.ξ) ** mat.C * (
        ramp(τ / mat.τ_base - mat.τ_cut) ** mat.m
    )

    ΔF_v = γ_dot * inv(F_e) * dev_σ_B / τ * F_e * F_v * dt
    return ΔF_v


##### defining mesh and boundary facets ############
l_x, l_y, l_z = 2, 4, 20

msh = mesh.create_box(
    MPI.COMM_WORLD,
    [[0.0, 0.0, 0.0], [l_x, l_y, l_z]],
    [2 * l_x, 2 * l_y, 2 * l_z],
    mesh.CellType.tetrahedron,
)


def bottom(x):
    return np.isclose(x[2], 0)


def top(x):
    return np.isclose(x[2], l_z)


tdim = msh.topology.dim  # topological dimension
fdim = msh.topology.dim - 1  # topological dimension of facets

bottom_facets = mesh.locate_entities_boundary(msh, fdim, bottom)
top_facets = mesh.locate_entities_boundary(msh, fdim, top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
    msh, 2, marked_facets[sorted_facets], marked_values[sorted_facets]
)


# defining function spaces
V = VectorFunctionSpace(msh, ("CG", 1))
T = TensorFunctionSpace(msh, ("CG", 1))


def interpolate_expression(tensor):
    expr = Expression(tensor, T.element.interpolation_points())
    interpolated = Function(T)
    interpolated.interpolate(expr)
    return interpolated


# defining measures
ds = Measure(
    "ds",
    domain=msh,
    subdomain_data=facet_tag,
    metadata={"quadrature_degree": 4},
)

dx = Measure(
    "dx",
    domain=msh,
    metadata={"quadrature_degree": 4},
)


########### boundary conditions ###############
u_bc1 = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)
bottom_dofs = locate_dofs_topological(V, fdim, facet_tag.find(1))
bc1 = dirichletbc(u_bc1, bottom_dofs, V)

u_bc2 = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)
top_dofs = locate_dofs_topological(V, fdim, facet_tag.find(2))
bc2 = dirichletbc(u_bc2, top_dofs, V)

bcs = [
    bc1,
    bc2,
]

tdim = msh.topology.dim  # topological dimension
fdim = msh.topology.dim - 1  # topological dimension of facets

δ = Identity(tdim)

def compute(mat: Material_BB, filename: str, iterations: int):
    file = XDMFFile(MPI.COMM_WORLD, filename, "w")
    file.write_mesh(msh)

    δu = TestFunction(V)
    du = TrialFunction(V)
    u = Function(V)
    u.name = "u"
    u0 = Function(V)
    u00 = Function(V)

    F = grad(u) + δ  # deformation gradient tensor
    J = det(F)  # Jacobian of F
    F_e = F
    F_v = F
    ΔF_v = as_tensor(
        ScalarType(
            [
                [0.0, 0.0, 0.0],
                [0.0, 0.0, 0.0],
                [0.0, 0.0, 0.0],
            ]
        )
    )
    time_values = np.linspace(0, iterations - 1, iterations)
    dt = time_values[1] - time_values[0]

    traction = Constant(msh, ScalarType([0.0, 0.0, 0.0]))
    f = Constant(msh, ScalarType([0.0, 0.0, 0.0]))

    σ_A = get_sigma_A(F, mat)
    σ_B = get_sigma_B(F_e, F_v, mat)

    σ = σ_A + σ_B

    i, j, k = indices(3)

    Form = mat.ρ * (u[i] - 2 * u0[i] + u00[i]) / dt**2 * δu[i] * dx
    Form += F[i, k] * σ[j, k] * δu[i].dx(j) * dx
    Form -= mat.ρ * f[i] * δu[i] * dx
    Form -= traction[i] * δu[i] * ds(2)

    Gain = derivative(Form, u, du)

    problem = NonlinearProblem(Form, u, bcs, Gain)
    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-4

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "preonly"
    opts[f"{option_prefix}pc_type"] = "lu"
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
    ksp.setFromOptions()

    #stress_values = []
    start = time()
    for n, t in enumerate(time_values):
        print(f"Iteration: {n+1}")
        print(f"Time of the last loop: {(time() - start):.2f}")
        start = time()
        F = grad(u) + δ
        J = det(F)
        F_v += ΔF_v
        F_e = F * inv(F_v)

        σ_A = get_sigma_A(F, mat)
        σ_B = get_sigma_B(F_e, F_v, mat)

        σ = σ_A + σ_B
        # traction.value = ScalarType([0.0, 0.0, 100*t])

        u_bc2 = np.array([0.0, 0.0, 0.001 * t * l_z], dtype=PETSc.ScalarType)

        problem.bcs = [
            bc1,
            dirichletbc(u_bc2, top_dofs, V),
        ]
        solver.solve(u)
        u.x.scatter_forward()

        interpolate_expression(σ)
        # stress_values.append(stress.x.array)

        ΔF_v = get_delta_F_v(F_v, F_e, σ_B, mat, dt)

        u00.x.array[:] = u0.x.array
        u0.x.array[:] = u.x.array
        file.write_function(u, n)

    file.close()
    #return stress_values


mat = Material_BB(
    2.05862,  # μ
    9.22561,  # λ_lock
    76.269,  # κ
    6.63277,  # s
    7.56263e-6,  # ξ
    -1.99759,  # C
    9.84628,  # τ_base
    22.3204,  # m
    0.01,  # τ_cut
    1,  # ρ
)

sigma_values = compute(mat, "test.xdmf", 100)

The output on the terminal:

Iteration: 61
Time of the last loop: 0.68
Iteration: 62
Time of the last loop: 0.65
Iteration: 63
Time of the last loop: 0.52
Iteration: 64
Time of the last loop: 0.54
Iteration: 65
Time of the last loop: 59.11
Iteration: 66
Time of the last loop: 61.88
Iteration: 67
Time of the last loop: 61.55
Iteration: 68
Time of the last loop: 64.58
Iteration: 69
Time of the last loop: 74.01

Apologies again for not sticking to the rules.

2

The issue in this code is that the expression that you are trying to interpolate gets more and more complicated over time (F_v is a sum that grows in time-dependent on itself and a highly nonlinear operator, get_delta_F_v).

You should aim for not having to compute a dolfinx.fem.Expression inside a temporal loop, as this causes JIT compilation of generated code.

I do not have any experience with the mathematical model that you are using. Maybe some else has the time to parse your code or has tried to do something similar. I would suggest trying to reduce the complexity of your problem, as most people do not have time to look through long codes trying to backward engineer what you have implemented.

2 Likes
3

You’ve helped me enormously as it is, I just needed somebody to point me in the right direction, thank you so so much.

This forces me think about the problem in a completely new way!

The problem is as follows - Bergstrom-Boyce derived a material model to describe elastomeric polymers (rubber etc.), in which they formulated the relationships between the stress tensor \sigma, the deformation gradient tensor F and the derivative of the deformation gradient tensor \dot{F}.

You can also derive from the balance of momentum the following expression:

\rho \cdot \frac{u^n- 2u^{n-1}+u^{n-2}}{\Delta t^{2}} - inv(F)\cdot\sigma\cdot det(F)-\rho\cdot f=0

which is always correct assuming F and \sigma are appropriate. I used the expression as my form and updated F and \sigma foe each time step as they change in time. Mathematically it checks out, simply not being familiar with the innards of Fenics, I didn’t think it would pose problems. I wonder why the solver has no trouble solving for u if the expression grows larger and larger at each iteration.

Ah well, back to the drawing board. Thanks a lot!

4

As far as I can tell you do not use sigma in your residual Form. Thus at each time step (n) your are solving Form(u_n)=0, and since you are solving a non-linear problem, the initial guess for the Newton iteration is u_{n-1}.

This means that the form does not grow over time.

I think that in general your sigma isn’t correct, as
Pre loop
Sigma=sigma_a(F(u_0))+sigma_b(F(u_0),F(u_0))

First step in loop
Sigma = sigma_a(F(u_1))+sigma_b(F(u_1)inv(F(u_1)), F(u_1))

Then things get complicated, as you add in delta_F_v which is always a function of u_i for any step in the loop after the first.

I am not sure if this is what you mean to do, or if sigma should be a tensor which you keep adding a contribution to for each time step i

1 Like
5

F(u_1) is still the identity tensor at this point, so it shouldn’t matter?

But thank you for going deeper into my code. I need to reassess my assumptions how Fenics calculates things under the hood.

6

In step 1 yes, But as I said, from step 2 an onwards this gets really complicated, where you add more and more terms (that all depend on the current time step only, not the previous loading conditions).

1 Like
7

My formulation of the problem stemmed from the fact that I didn’t understand how Fenics handles expression objects and how they’re summed and injected into the form etc.

But the idea was that one of the constitutive equations was.
\dot{F_v^n} = \gamma\cdot (F_e^n)^{-1}AF_e^nF_v^n
which is roughly equal to:
\frac{\Delta{F_v^n}}{\Delta t} = \gamma\cdot (F_e^n)^{-1}AF_e^nF_v^n
which gives
\Delta F_v^n = \gamma\cdot (F_e^n)^{-1}AF_e^nF_v^n\cdot\Delta t
which finally gives:
F_v^{n+1} = F_v^n+\Delta F_v^n
so that was the logic behind updating those expressions.

Thanks again for taking the time to look at my problem and giving some feedback. :pray: