Cahn-Hilliard equation with a function as the initial guess of the Newton solver

Question

Hi folks,

Just a quick question about Cahn-Hilliard equation. Consider (from demo):

# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self):
        random.seed(2 + MPI.process_number())
    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,)

Let us assume that instead of "0.63 + 0.02*(0.5 - random.random())" for values[0], we want to assign a function instead. In other words, assume that our initial guess for the Newton solver is given by a function.

I would be grateful if you could let me know how to implement it.

Cheers,
Mo

Answer 1

Hi,

If you for example have a function of the coordinates, func, then you could do something like this

def func(x):
   return x

class InitialConditions(Expression):
  def eval(self, values, x):
    values[:] = func(x)
  def value_shape(self):
    return (2,)

Comments

Hi Mikael,

Thank you for your reply. My function is a vector valued function. More clearly, if our mixed function space is ME=V*V (similar to the demo example), my vector valued function is the solution of an ODE defined on V.

Sorry for not being clear at first.

Cheers,
Mo

---

Do you mean that you have a dolfin Function on V that you want to use for initializing your solution Function on ME? Then you don't need an Expression. You simply want to assign the V Function to one of the two components in ME?

---

That's exactly what I need Mikael. Would you please let me know how to implement that? Let us consider the following simplified demo.

import random
from dolfin import *

# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self):
        random.seed(2 + MPI.process_number())
    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
        self.reset_sparsity = True
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity)
        self.reset_sparsity = False

# Create mesh and define function spaces
mesh = UnitSquareMesh(96, 96)
V = FunctionSpace(mesh, "Lagrange", 1)
ME = V*V

# 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()
u.interpolate(u_init)
u0.interpolate(u_init)

u_t = Function(V)
u_t.vector()[:] = 1

u.sub(0) = u_t

c,  mu  = split(u)

c_t = project(c,V)

print c_t.vector().array()  

Let's say we want to assign u_t to c. Obviously, u.sub(0) = u_t does not work and I get the following error:

  File "CH.py", line 53
    u.sub(0) = u_t
SyntaxError: can't assign to function call

Thank you again for helping me out.

Cheers,
Mo

---

Ok,
If you have a development of dolfin, then you can use the FunctionAssigner class. But you probably have an older version? You can then assign to the solution in ME using the dofmap, like this to assign to the first subspace of ME (from memory, not tested)

u = Function(ME)
u1 = Function(V)
vals = ME.sub(0).dofmap().collapse(mesh)[1].values() # dofs of first subspace
u.vector()[vals] = u1.vector()[:]

---

Thank you so much Mikael, this works fine.

Cheers,
Mo