No convergence for outsourced function which is included in the variation formulation

hunterr · May 31, 2023, 12:09pm

Hello everyone,

i am currently writing my first python package for modeling a hyperelastic problem. Therefor I want to write the different material laws and kinematic quantities in separate files and call them from a central file.
Unfortunately, the solution does not converge if the ufl variables for the kinematic quantities are not defined directly at the weak form of the problem itself.

Since the problem cannot be easily translated into a short MWE, I’ll try the following example.
I assume that I have fundamentally misunderstood or misapplied something here and the error can be easily detected by someone more experienced than me.

Code that works:

# FILE A: material_models
import ufl
from dolfinx import fem
from petsc4py.PETSc import ScalarType

def neo_hooke(domain, u):
    I = ufl.variable(ufl.Identity(len(u)))
    F = ufl.variable(I + ufl.grad(u))
    C = ufl.variable(F.T * F)
    I_1 = ufl.variable(ufl.tr(C))
    J = ufl.variable(ufl.det(F))
    E = ScalarType(1.0e4)
    nu = ScalarType(0.3)
    mu = fem.Constant(domain, E/(2*(1 + nu)))
    lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
    psi = (mu / 2) * (I_1 - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
    P = ufl.diff(psi, F)
    S = 2 * ufl.diff(psi, C)
    return P, S

Code that does not work:

# FILE A: material_models

import ufl
from dolfinx import fem
from petsc4py.PETSc import ScalarType
from kinematics import *

def neo_hooke(domain, u):
    E = ScalarType(1.0e4)
    nu = ScalarType(0.3)
    mu = fem.Constant(domain, E/(2*(1 + nu)))
    lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
    psi = (mu / 2) * (I1(u) - 3) - mu * ufl.ln(jacobian_defgrad(u)) + (lmbda / 2) * (ufl.ln(jacobian_defgrad(u)))**2
    P = ufl.diff(psi, deformation_gradient(u))
    S = 2 * ufl.diff(psi, right_cauchygreen(u))
    return P, S
    
# FILE B: kinematics
import ufl
import numpy as np

def identity_tensor(u):
    d = len(u)
    return ufl.variable(ufl.Identity(d))

def deformation_gradient(u):
    I = identity_tensor(u)
    return ufl.variable(I + ufl.grad(u))

def jacobian_defgrad(u):
    F = deformation_gradient(u)
    return ufl.variable(ufl.det(F))

def right_cauchygreen(u):
    F = deformation_gradient(u)
    return ufl.variable(F.T * F)

def I1(u):
    C = right_cauchygreen(u)
    I1 = ufl.tr(C)
    return ufl.variable(I1)

Thanks in advance!

dokken · May 31, 2023, 1:19pm

Why cant you make this into an MWE?

Simply create a unit square or cube
Define the function space V and the function u.
Create the variational form in 2 ways.

This follows here:

import ufl
from dolfinx import fem, mesh
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
def neo_hooke(domain, u):
    I = ufl.variable(ufl.Identity(len(u)))
    F = ufl.variable(I + ufl.grad(u))
    C = ufl.variable(F.T * F)
    I_1 = ufl.variable(ufl.tr(C))
    J = ufl.variable(ufl.det(F))
    E = ScalarType(1.0e4)
    nu = ScalarType(0.3)
    mu = fem.Constant(domain, E/(2*(1 + nu)))
    lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
    psi = (mu / 2) * (I_1 - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
    P = ufl.diff(psi, F)
    S = 2 * ufl.diff(psi, C)
    return P, S


def neo_hooke2(domain, u):
    E = ScalarType(1.0e4)
    nu = ScalarType(0.3)
    mu = fem.Constant(domain, E/(2*(1 + nu)))
    lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
    psi = (mu / 2) * (I1(u) - 3) - mu * ufl.ln(jacobian_defgrad(u)) + (lmbda / 2) * (ufl.ln(jacobian_defgrad(u)))**2
    P = ufl.diff(psi, deformation_gradient(u))
    S = 2 * ufl.diff(psi, right_cauchygreen(u))
    return P, S
    
def identity_tensor(u):
    d = len(u)
    return ufl.variable(ufl.Identity(d))

def deformation_gradient(u):
    I = identity_tensor(u)
    return ufl.variable(I + ufl.grad(u))

def jacobian_defgrad(u):
    F = deformation_gradient(u)
    return ufl.variable(ufl.det(F))

def right_cauchygreen(u):
    F = deformation_gradient(u)
    return ufl.variable(F.T * F)

def I1(u):
    C = right_cauchygreen(u)
    I1 = ufl.tr(C)
    return ufl.variable(I1)


domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.VectorFunctionSpace(domain, ("Lagrange", 1))
u = fem.Function(V)


neo_hooke(domain, u)

neo_hooke2(domain, u)

Neither of these produce an error for me, so there is no way of debugging your issue.

hunterr · May 31, 2023, 2:36pm

MWE here we go:

FILE 1: test.py

import numpy as np
import ufl
from dolfinx import fem, mesh, log
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from dolfinx import nls
from kinematix import *

# domain geometries
length = 0.25
height = 0.02

# dimensions of indenter 
D = 0.15
indentation_depth = height/2.0
center_x = length / 2
center_y = (height + (D/2)) - indentation_depth

domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0],[length, height]], [40, 9])
x = ufl.SpatialCoordinate(domain)
n = ufl.FacetNormal(domain)

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

def top(x):
    return np.isclose(x[1], height)

tdim = domain.topology.dim
fdim = tdim - 1

bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
top_facets = mesh.locate_entities_boundary(domain,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(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

metadata = {"quadrature_degree":4}
ds = ufl.Measure('ds',domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
V = fem.VectorFunctionSpace(domain, ("CG",1))

u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_tag.find(1)), V)

# definition of gap function 
def gapp(u, x):
    R = ufl.sqrt(pow(((x[0] - center_x) + u[0]), 2) + pow(((x[1] - center_y) + u[1]), 2)) 
    return  R - (D/2)

# maculay bracket
def maculay(x): 
    return (x+abs(x))/2

# kinematic quantities
I = ufl.variable(ufl.Identity(len(u)))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
I_1 = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Neo-Hooke
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))

################################## questionable part ##################################
psi = (mu / 2) * (I_1 - 3)   - mu * ufl.ln(J)                   + (lmbda / 2) * (ufl.ln(J))**2 # works
#psi = (mu / 2) * (I1(u) - 3) - mu * ufl.ln(jacobian_defgrad(u)) + (lmbda / 2) * (ufl.ln(jacobian_defgrad(u)))**2 # does not work
#######################################################################################

# 1st Piola-Kirchhoff
P = ufl.diff(psi, F)

# Penalty Stiffness
kpen = fem.Constant(domain, ScalarType(1e7))

# penalty term for contact
penalty = kpen * n *maculay(-gapp(x, u))

F1 = ufl.inner(ufl.grad(u_), P)*dx
F2 = kpen * ufl.inner(u_, n*maculay(-gapp(x, u)))*ds(2)

form = F1 + F2
Jac = ufl.derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=Jac)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

solver.atol = 1e-6
solver.rtol = 1e-6
solver.convergence_criterion = "incremental"
solver.report = True
ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)

w, converged = solver.solve(u)
assert(converged)
print(f"Number of interations: {w:d}")

with XDMFFile(domain.comm, "mwe_package/result/output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u)

FILE 2: kinematix.py

import ufl

def identity_tensor(u):
    d = len(u)
    return ufl.variable(ufl.Identity(d))

def deformation_gradient(u):
    I = identity_tensor(u)
    return ufl.variable(I + ufl.grad(u))

def jacobian_defgrad(u):
    F = deformation_gradient(u)
    return ufl.variable(ufl.det(F))

def right_cauchygreen(u):
    F = deformation_gradient(u)
    return ufl.variable(F.T * F)

def I1(u):
    C = right_cauchygreen(u)
    I1 = ufl.tr(C)
    return ufl.variable(I1)

hunterr · June 2, 2023, 5:52am

Is it possible for you to reproduce this issue @dokken ?

Topic		Replies	Views
Define Nonlinear problem with ufl.derivative and MixedElement variational formulation dolfinx	1	42	April 18, 2024
How to define a constant in variational formulation? General	1	53	April 9, 2024
Calculate adjoint derivatives w.r.t to multiple Constants dolfinx	6	335	December 4, 2023
A demo from MHD General	0	225	October 10, 2023
Convergence in Stokes problem variational formulation	2	229	August 8, 2023

No convergence for outsourced function which is included in the variation formulation

Related Topics