Modifying Cahn-Hilliard tutorial by adding Fischer-Burmeister function

dolfinx
1

Hello everyone,

I have been trying to add Fischer-Burmeister functions (Nonlinear Complementary functions) to the existing Cahn-Hilliard tutorial problem to impose a constrain on the concentration value p to be within range of 0 to 1. It is needed to be solved along with the set of equation but the constrain functions were not in any integral. Hence directly adding these along with the tutorial equation is making error as


Traceback (most recent call last):

File "/mnt/e/PPP/Project/Chapter 1-29-11/Phase_field solver/code/Cahn_hilliard_h.py", line 57, in <module>

F = F0 + F1 + F2 + F3

TypeError: unsupported operand type(s) for +: 'Form' and 'Sum'

The set of equations to be solved were look like


F0 = inner(p, q) * dx - inner(p0, q) * dx + dt * k * inner(grad(mu_mid), grad(q)) * dx

F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - 3*epsilon*areainteng*inner(grad(p), grad(v)) * dx + inner(lam1,v)*dx - inner(lam2,v)*dx

F2 = (ufl.sqrt((p-1)**2+lam1**2) + (p-1) - lam1)

F3 = (ufl.sqrt((-p)**2+lam2**2) - p - lam2)

Additionally i also did changes to the functions as i needed to find the value of lamda (lam) for all the node including with the chemical potential and concentration. The modification is as follows


#Mesh creation

msh = create_unit_square(MPI.COMM_WORLD,24,24, CellType.quadrilateral)

P1 = element("Lagrange", msh.basix_cell(), 1)

ME = functionspace(msh, mixed_element([P1, P1, P1, P1]))

#Test functions

q, v, l1, l2 = ufl.TestFunctions(ME) #q-chemical potential,v-order parameter

u = Function(ME) # current solution

u0 = Function(ME) # solution from previous converged step

# Split mixed functions

p, mu,lam1,lam2= ufl.split(u)

p0, mu0,lam01,lam02= ufl.split(u0)

Is there any way to solve these set of equation. Also i have tried multiplying dx to these but the solver is making error. Please suggest me some approaches for handling the same. Any advice would be truly appreciated. Thanks for your time and consideration.

2

You need to copy the entire code, otherwise it’s impossible to understand where the error is just looking at those two snippets.

3

Sorry for the inconvenience. The entire code is


import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
import dolfinx

# Save all logging to file
log.set_output_file("log.txt")

lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
msh = create_unit_square(MPI.COMM_WORLD, 50, 50, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1, P1, P1]))
q, v, l1, l2 = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step

# Split mixed functions
c, mu, lam1, lam2 = ufl.split(u)
c0, mu0, lam01, lam02 = ufl.split(u0)

# Zero u
u.x.array[:] = 0.0

# Interpolate initial condition
np.random.seed(30)
u.sub(0).interpolate(lambda x: 0.75 + 0.02 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()

# Compute the chemical potential df/dc
c = ufl.variable(c)
f = 100 * c**2 * (1 - c)**2
dfdc = ufl.diff(f, c)

# mu_(n+theta)
mu_mid = (1.0 - theta) * mu0 + theta * mu

# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx

#NCP functions
F2 = (ufl.sqrt((c-1)**2+lam1**2) + (c-1) - lam1)
F3 = (ufl.sqrt((-c)**2+lam2**2) - c - lam2)

F = F0 + F1 + F2 + F3  #Adding for solving together the weak form with the functions

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver used inside the NewtonSolver by
# 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"
ksp.setFromOptions()

# Output file
file = XDMFFile(MPI.COMM_WORLD, "demotest.xdmf", "w")
file.write_mesh(msh)

# Step in time
t = 0.0
timesteps = [0]

# Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 3 * dt
else:
T = 50* dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector

V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array

while (t < T):
t += dt
r = solver.solve(u)
print(f"Step {int(t/dt)}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
print(u.x.array)

file.write_function(c, t)
file.close()

Here the equation to be solved with the two NCP functions (here F2 and F3) for all the elements as it acts as a constraint to make the value of c within range of 0 and 1 ever after n iterations. After considerable iterations the value of c is going negative, in order to avoid that, this needed to be solved.

Thanks again for your time and considerations.

4

You should try and make it as easy as possible for users of this forum to help you. The code that you posted is not indented correctly, and thus people will not be able to run it.

5

I am sorry for that. While pasted here, got alignment changed last time. I hope now it will work. Thanks for your support.

import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
import dolfinx

# Save all logging to file
log.set_output_file("log.txt")
lmbda = 1.0e-02  # surface parameter
dt = 5.0e-06  # time step
theta = 0.5  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
msh = create_unit_square(MPI.COMM_WORLD, 24, 24, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1, P1, P1]))
q, v, l1, l2 = ufl.TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step
# Split mixed functions
c, mu, lam1, lam2 = ufl.split(u)
c0, mu0, lam01, lam02 = ufl.split(u0)
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition
np.random.seed(30)
u.sub(0).interpolate(lambda x: 0.75  + 0.02 * (0.5 - np.random.rand(x.shape[1])))
#print (u.x.array)
u.x.scatter_forward()
# Compute the chemical potential df/dc
c = ufl.variable(c)
#f = 100 * c**2 * (1 - c)**2
f = 100 * c**2 * (1 - c)**2+((3/2)*(grad (c))**2)
dfdc = ufl.diff(f, c)
# mu_(n+theta)  
mu_mid = (1.0 - theta) * mu0 + theta * mu

# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F2 = (ufl.sqrt((c-1)**2+lam1**2) + (c-1) - lam1)
F3 = (ufl.sqrt((-c)**2+lam2**2) - c - lam2)
F = F0 + F1 + F2 + F3 

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver used inside the NewtonSolver by
# 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"
ksp.setFromOptions()

# Output file
file = XDMFFile(MPI.COMM_WORLD, "demotest.xdmf", "w")
file.write_mesh(msh)

# Step in time
t = 0.0
timesteps = [0]
#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
    T = 3 * dt
else:
    T =  150* dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array
while (t < T):
    t += dt
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
    print(u.x.array)
    file.write_function(c, t)

file.close()
6

F2 and F3 are certainly missing a test function. You haven’t written the actual weak formulation, so below I’ve taken a guess.

diff --git a/run.py b/run.py
index 7776c59..68af6b6 100644
--- a/run.py
+++ b/run.py
@@ -45,8 +45,8 @@ mu_mid = (1.0 - theta) * mu0 + theta * mu
 # Weak statement of the equations
 F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
 F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
-F2 = (ufl.sqrt((c-1)**2+lam1**2) + (c-1) - lam1)
-F3 = (ufl.sqrt((-c)**2+lam2**2) - c - lam2)
+F2 = (ufl.sqrt((c-1)**2+lam1**2) + (c-1) - lam1) * l1 * dx
+F3 = (ufl.sqrt((-c)**2+lam2**2) - c - lam2) * l2 * dx
 F = F0 + F1 + F2 + F3
 
 # Create nonlinear problem and Newton solver
@@ -89,4 +89,3 @@ while (t < T):
     print(u.x.array)
     file.write_function(c, t)
 
-
1 Like
8

Thanks for all your support. I have a small question, could anyone confirm that the implementation of the Nonlinear complementary conditions (Fischer Burmeister conditions)functions used as F2 and F3 in the above code is correct? I am little unsure whether implementation of the above necessarily incorporate as a constraint to the variable c. It will be very much helpful and thanks again for your time and understanding.

9

Not sure how many people in the discourse community are familiar with Fischer Burmeister conditions to be able to check. I certainly I am not :wink:

10

Thanks for your response. May I know, is there any approach in dolfinx that I can have my functions F2 and F3 to be applied to all the nodal points instead of ‘dx’ used here to specify the complete domain?. As I need to apply the function F2 and F3 to each of the nodal points only and not needed to be interpolated like in F0 and F1. Is there any specific way to describe my equations like that?
Or maybe I am wrong in my understanding about dx ?

Thanks for your time.