Getting error in Newton solver by adding elasticity terms in Cahn hilliard equation

General
,
1

Hello all, I am trying to incorporate elasticity terms in the chemical potential equation. I have added term3 and term4 in code provided below in the Cahn Hilliard expression. If I remove term3 and term4, solver is working perfectly. I am getting the error meesage “Applying nonlinear operator Power to expression depending on form argument v_1”.

Can someone please help me to figure out the solver issue, I guess I have to modify the Non-linear Newton solver to accomodate elasticity terms.

Code

comm = MPI.COMM_WORLD

msh = create_unit_square(MPI.COMM_WORLD, 96,96 , CellType.triangle)



# finite element function space
element_u = basix.ufl.element("Lagrange", msh.basix_cell(), degree=1, shape=(msh.geometry.dim,))
V = fem.functionspace(msh, element_u)

# Define the state
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Domain measure.
dx = ufl.Measure("dx", domain=msh)

# cahn hilliard function space

lmbda = 1.0e-02  # surface parameter
dt = 5.0e-08  # time step
theta = 1  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson

P1 = basix.ufl.elementelement("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
ME = fem.functionspace(msh, mixed_element([P1, P1]))

q1, v1 = ufl.TestFunctions(ME)

cu = Function(ME)  # current solution
cu0 = Function(ME)  # solution from previous converged step

# Split mixed functions
c, cmu = ufl.split(cu)
c0, cmu0 = ufl.split(cu0)

# initialize your concentrations
# Zero u
cu.x.array[:] = 0.0

# Interpolate initial condition
rng = np.random.default_rng(42)
cu.sub(0).interpolate(lambda x: 0.63 + 0.02 * (0.5 - rng.random(x.shape[1])))
cu.x.scatter_forward()

c = ufl.variable(c)
f = 100 * c**2 * (1 - c) ** 2
dfdc = ufl.diff(f, c)

# mu_(n+theta)
mu_mid = (1.0 - theta) * cmu0 + theta * cmu


# define the constants for elasticity

Ep = 3.3
Epo = 3.1

E = c*Epo + (1-c)*Ep
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
# Plane strain definition
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

# Identity tensor
I = ufl.Identity(2)


# Weak form components
F0 = (ufl.inner(c, q1) - ufl.inner(c0, q1) + 
      dt * ufl.inner(ufl.grad(mu_mid), ufl.grad(q1))) * ufl.dx

# Nonlinear terms
eps_u = ufl.sym(ufl.grad(u))  # Assuming u is defined elsewhere
trace_eps = ufl.tr(eps_u)
term3 = inner((nu * (Epo - Ep) * (trace_eps**2) / (2 * (1 + nu) * (1 - 2 * nu))), v1) *dx
term4 = inner(((Epo - Ep) / (2 * (1 + nu))) * ufl.tr(eps_u * eps_u), v1) * dx

# Weak form F1
F1 = inner(mu, v1) * dx - inner(dfdc, v1) * dx - lmbda * inner(grad(c), grad(v1)) * dx - term3 - term4

# Total weak form
F2 = F0 + F1


# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F2, cu)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-5
solver.max_it = 200 # Maximum number of iterations
# We can customize the linear solver used inside the NewtonSolver 

# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
sys = PETSc.Sys()  # type: ignore
# For factorisation prefer superlu_dist, then MUMPS, then default
if sys.hasExternalPackage("superlu_dist"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
elif sys.hasExternalPackage("mumps"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
2

You can’t use TrialFunction when you want to solve a non-linear problem.

This variable should be your solution function cu.

3

What I’m attempting to convey here is that the displacement (u) is dependent on the concentration component (c) of cu. In my case their is a correlation between the strain and the concentration. So, I will solve the elasticity equation and cahn hilliard equation in the staggered approach. Maybe, I have to define the same function space for both the equation?

4

I replace the u with cu and I didn’t get any error and it’s not right as well. But my point is that as concentration (c) is dependent on u (displacement), so how to proceed with that.

5

If you have two unknowns, concentration and displacement, that you need to solve for, you need to create a mixed problem. See for instance the Stokes demo in the dolfinx documentation, or something closer to your problem as it is non-linear:
http://jsdokken.com/FEniCS-workshop/src/multiphysics/coupling.html

6

ok, so what I am doing is from the Cahn hilliard equation solver which is available in fenics demo section, i am solving for concentration (c) and chemical potential (cmu). Then, from the elasticity solver available, I am solving for displacements. So, total I have 3 unknowns: - concentration, chemical potential and displacements. I am solving in staggered approach. I guess if I define mixed element function in which concentration and chemical potential are scaler quantities from Cahn hilliard equation and displacement which is vector quantity of dimension 2 from elasticity equation, will the solver works?

Thanks for helping and understanding my problem.

7

If you are solving in a staggered approach, then
You should incorporate the solution of the displacement in your variational form, not the total function of that space.

8

Hi, thanks for the help. So, I define u as a function instead of a Trail function and it worked. I read that we can’t use Trial function in a non-linear form.
But, my solution is not-converging. I increased my maximum iterations to 1000 but still its not converging. So, I am guessing, the way I have implemented the staggered approach is wrong. Any suggestions?

9

Please make an updated minimal example that illustrates your approach, and point to which of the equations that does not converge.

10

sure, I am working on the code to see any bugs, thanks.