Define an expression based on current solution

variational formulation
1

What is the best way to define an Expression (in c code) which makes use of the current solution?

Working from the Cahn-Hilliard demo where the function u is split into ‘c’ and ‘mu’, the link suggests:
dfdc = Expression('200*c*(2*pow(c,2)-3*c+1)',c = u.sub(0), degree=1)

which works, but appears slow when implemented in a modified CH demo below. Note I also implement the Jacobian contribution as an expression. Flag = True uses expressions for dfdc and its derivative, Flag = False is the original form. I also noted more Newton steps which makes me wonder if it is related to the Jacobian contribution?
Any help understanding this / improving it would be appreciated. Thanks.

import random
from dolfin import *

flag = False


# Class representing the intial conditions
class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)

# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

# Model parameters
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-Nicolson

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

# Create mesh and build function space
mesh = UnitSquareMesh.create(96, 96, CellType.Type.triangle)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)
# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)
if flag:
    dfdc = Expression('200*c*(2*pow(c,2)-3*c+1)',c = u.sub(0), degree=1)
    dfdcd = Expression('200*(1-6*c+6*pow(c,2))',c = u.sub(0), degree=1)
dfdcv = variable(dfdc)

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

# Weak statement of the equations
L0 = c*q - c0*q + dt*dot(grad(mu_mid), grad(q))
L1 = mu*v - dfdcv*v - lmbda*dot(grad(c), grad(v))
F = L0 + L1
L = F*dx

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# derivative of dfdc has no information as expression.
if flag:
    a += dc*dfdcd*(diff(F,dfdcv))*dx

# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("output.pvd", "compressed")

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    file << (u.split()[0], t)
1 Like
2

Update: I increased the ‘degree=2’ in the expression definition and that fixed the poor newton method convergence, but it is still significantly slower.

3

Hi Mike, I have been facing similar issues with convergence of the Newton method especially when using flags (‘conditionals’ in my case). I believe that it is related to the Jacobian as you have mentioned, since it needs to be ‘exact’ in order for the Newton’s method to achieve quadratic convergence. So, when you increase the order of the expression, the accuracy of the interpolation onto each element in the space is increased. This should be the reason for the improved convergence. But you say that it is still significantly slower. Compared to what? Are you comparing the two cases when flag=true and flag=false ? Or when there is no flag at all?
I am working in python. So, I am using a ‘conditional’ function to switch between definitions of a give function. For some reason boolean statements like ‘if’ donot work when evaluating UFL form. In my case the ‘conditional’ should switch between whether the PDE being solved is linear or nonlinear and in any case I solve with the Nonlinear solver. And I compare this with when I have no ‘conditional’ at all and I directly solve the nonlinear equations. The second case is much faster even when the first case is solving the same nonlinear equations. What I want to say is, the evaluation of the ‘conditional’ itself is taking a lot of time. Maybe you can check the same by removing the flag altogether and try solving the both cases separately.

For your original question on what is the best way to define an expression using current solution, I am still looking for an answer. Please let me know if come across something even if it is in c code.

Thanks,
Sid