Newton Solver for Cahn-Hilliard-Stokes problem not converging

Linear Algebra
1

Dear All,
I’m currently trying to solver the coupled system of Cahn-Hilliard alongside the Stokes equation.

Screen Shot 2022-10-20 at 17.56.23
Screen Shot 2022-10-20 at 17.56.231274×484 27.3 KB

I tried to solve it using a mixed function space method where is have elements for [conc,chemical_potenital,pressure,velocity] , and as a first test i made it completely decoupled by not having the coupling terms. However somewhere along the 4th time step of the newton solver intially converges and then diverges. I do not understand the source of the problem , as i am right now just solving two decoupled problems using mixed formulation. I have attached the code below . Would be gratefull to recieve any form of help :slight_smile:

import dolfin as df
import matplotlib.pyplot as plt
import numpy as np
from ufl import tanh
import sys

a=1
b=1e-2
l_cap=np.sqrt(a)
M=1
dt=5.0e-04
x=0.5
num_steps=10
nt=64

mesh=df.UnitSquareMesh(nt,nt)

boundary_markers = df.MeshFunction("size_t",mesh, mesh.topology().dim()-1,0)
regions=df.MeshFunction("size_t",mesh, mesh.topology().dim())

class BoundaryBottom(df.SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and df.near(x[1], 0, self.tol)



class BoundaryTop(df.SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and df.near(x[1], 1, self.tol)



class PeriodicBoundary(df.SubDomain):

    def inside(self, x, on_boundary):
        return bool (df.near(x[0], 0) and on_boundary)

    def map(self, x, y):
        y[0]=x[0]-1
        y[1]=x[1]



def CH_form():
    c_form= c_1*phi -c_0*phi + dt*M*df.inner(df.grad(mu_0*(1-x)+mu_1*x),df.grad(phi))
    mu_form=mu_1*lam -(2*a*c_1*(1-c_1)*(1-2*c_1))*lam -b*df.inner(df.grad(c_1),df.grad(lam))
    form= (c_form +mu_form)*df.dx
    return form


def stokes_form():
	form=(df.inner(df.grad(u), df.grad(v)) + df.div(v)*p + q*df.div(u))*df.dx
	return form



Bottom=BoundaryBottom()
Top=BoundaryTop()
Left=BoundaryLeft()
Right=BoundaryRight()

pbc=PeriodicBoundary()


SFS = df.FunctionSpace(mesh, 'P', 1,constrained_domain=pbc)
s_1_element= df.FiniteElement('P', mesh.ufl_cell(),1)
s_2_element= df.FiniteElement('P', mesh.ufl_cell(),1)
p_element=df.FiniteElement('P', mesh.ufl_cell(),1)
v_element=df.VectorElement('P',mesh.ufl_cell(),2)

element = df.MixedElement([s_1_element, s_2_element,p_element,v_element])
FS=df.FunctionSpace(mesh,element,constrained_domain=pbc)

phi,lam,q, v = df.TestFunctions(FS)


f_1=df.Function(FS)
f_0=df.Function(FS)

c_0, mu_0,p_0,u_0=df.split(f_0)
c_1, mu_1,p,u=df.split(f_1)
#c_0=df.interpolate(df.Expression('1+tanh(((pow(((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5)),0.5)-0.2)/l))',l=l_cap,degree=1),FS.sub(0).collapse())
c_0=df.interpolate(df.Expression('0.5',degree=1),FS.sub(0).collapse())
noise =  0.1*(2*np.random.random(c_0.vector().size())-1)
noise -= noise.mean()
c_0.vector()[:] = c_0.vector()[:] + noise
df.plot(c_0)
plt.show()
#sys.exit()

f_0.assign(f_1)

form=CH_form() + stokes_form()

top_vel = df.Expression(("1", "0.0"),degree=1)
bottom_vel=df.Expression(("-1", "0.0"),degree=1)

bc1 = df.DirichletBC(FS.sub(3), top_vel, Top)
bc2=df.DirichletBC(FS.sub(3), bottom_vel, Bottom)
bcs=[bc1,bc2]

ConcFile=df.XDMFFile('conc.xdmf')
VelFile=df.XDMFFile('vel.xdmf')
t=0
for steps in range(num_steps):
    t+=dt
    df.solve(form==0,f_1,bcs,solver_parameters={'newton_solver': {'relaxation_parameter':1e0,'relative_tolerance':1e-7}})
    f_0.assign(f_1)
1 Like