Hi All,
Am using the Cahn-hilliard demo to solve a system of nonlinear pdes on mixed function space. But am not sure where to apply the Dirichlet boundary conditions. Also how should the form of my time stepping loop be? Am not sure of using the .vector() because my function space already has two vectors in it. I hope my question is clear enough.
A sample of my code is below.
The form of my functions: U, u, phi, p.
U and u are vector functions, and phi and p are scalar functions
Thank you to anyone that helps.
I have the boundary condition as list in bcs
from dolfin import *
class InitialConditions(Expression):
def eval(self, values, x):
values[0] = x[1]
values[1] = x[1]
values[2] = x[1]
values[3] = x[1]
values[4] = 0.01
values[5] = 0
def value_shape(self):
return (6,)
class MyEquations(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
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
mesh = RectangleMesh(Point(-1,-1),Point(1,1),200,200) #Domain[-1,1]*[-1,1]
V1= VectorFunctionSpace(mesh, 'CG', degree=2)
V2=VectorFunctionSpace(mesh, 'CG', degree=2)
Q1 = FunctionSpace(mesh, 'CG', degree=1 )
Q2 = FunctionSpace(mesh, 'CG', degree=1)
VQ = MixedFunctionSpace([V1,V2, Q1, Q2])
# Define trial and test functions
duu = TrialFunctions(VQ)
m,h,q,r = TestFunctions(VQ)
H= Function(VQ) # Current solution
H0 = Function(VQ) # solution from previous converged steps
dU, du, dphi, dp = split(duu)
U, u, phi, p = split(H)
U0, u0, phi0, p0 = split(H0)
#Create initial conditions and interpolate
u_init = InitialConditions()
H.interpolate(u_init)
H0.interpolate(u_init)
# Dirichlet boundary condition
# Sub domain for Dirichlet boundary condition
class Top(SubDomain):
def inside(self, x, on_boundary):
return (x[1]>1 - DOLFIN_EPS and on_boundary)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return (x[1]< -1 + DOLFIN_EPS and on_boundary)
UB= Constant((-1,0))
UT= Constant((1,0))
ub= Constant((-1,0))
ut= Constant((1,0))
top= Top()
bott = Bottom()
bcUB= DirichletBC(VQ.sub(0), UB ,bott)
bcUT= DirichletBC(VQ.sub(0), UB ,top)
bcub= DirichletBC(VQ.sub(1), ub ,bott )
bcut= DirichletBC(VQ.sub(1), ub ,top )
bcs = [bcUB,bcUT, bcub,bcut]
L = L0 + L1
# Is it ok to apply the boundary conditions here?
a = derivative(L,H, duu, bcs) # boundary conditions here?
problem = MyEquations(a, L)
solver = NewtonSolver("lu")
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
file = File("output.pvd", "compressed")
# Not very sure about the form here since Function space already has two vectors in it
# Is it right to use the .vector() here to copy H to H0 at the beginning of each time step?
# because I have mixed function space with two vector functions and two scalar functions
T_total = 50*dt
t=0.0
while (t<T_total):
t+=dt
H0.vector()[:] = H.vector() # Is it right to use the .vector( ) here ?
solver.solve(problem, H.vector())
