Picard Iteration Not Updating Properly

1

Hey everyone, I’m new here. I’m having an issue trying to solve the steady NS equation with a Picard Iteration. I’m trying to just solve the pressure driven flow through a channel with an incline. After one iteration, the solution appears to solve to the stokes flow result, which is expected from the zero initial guess being multiplied in the convective term, but the next iteration ends up solving to the same flow state and the error between the new and previous solution becomes zero.

Pretty new to coding up these types of problems, and I just can’t seem to see where I went wrong. Any help would be appreciated, thank you!

import fenics; import numpy; import sympy; import time; import dolfin; import mshr
import matplotlib.pyplot as plt

##########################################################################
# Define Geometry
Hf = 0.2; Hr = Hf/4;
L = 1; Lf = 0.25; Lr = 0.25;

point0 = fenics.Point(0,0)
point1 = fenics.Point(Lf,0)
point2 = fenics.Point(L-Lr,Hf-Hr)
point3 = fenics.Point(L,Hf-Hr)
point4 = fenics.Point(L,Hf)
point5 = fenics.Point(0,Hf)

domain = mshr.Polygon([point0, point1, point2, point3, point4, point5])

##########################################################################
# Generate Mesh
N_bulk = 50
mesh = mshr.generate_mesh(domain, N_bulk)

##########################################################################
# Fluid Properties
mu = 8.9e-4         # dynamic viscosity
rho = 1e3           # density
P_left = 10; P_right = 0

f   = fenics.Constant((0, 0))
mu  = fenics.Constant(mu)
rho = fenics.Constant(rho)
p_left = fenics.Constant(P_left)
p_right = fenics.Constant(P_right)

##########################################################################
# Define Function Spaces
# Velocity Elements
u_type = 'CG'; u_order = 2
u_element = fenics.VectorElement(u_type, mesh.ufl_cell(), u_order)
u_space = fenics.FunctionSpace(mesh, u_element)

# Pressure Elements
p_type = 'CG'; p_order = 1
p_element = fenics.FiniteElement(p_type, mesh.ufl_cell(), p_order)
p_space = fenics.FunctionSpace(mesh, p_element)

# Mixed Element and Function Space - (P2 / P1) Taylor Hood Elements
mixed_element = fenics.MixedElement([u_element, p_element])
mixed_space = fenics.FunctionSpace(mesh, mixed_element)

##########################################################################
# Define Trial Functions
u, p = fenics.TrialFunctions(mixed_space)

# Define Test Functions
v, q = fenics.TestFunctions(mixed_space)

# Define Previous Iteration
Wn = fenics.Function(mixed_space)
un, pn = Wn.split()

# Define Solution Function
W = fenics.Function(mixed_space)

##########################################################################
# Define Boundaries
# Inlet Surface
class Inlet(fenics.SubDomain):
	def inside(self, x, on_boundary):
		return fenics.near(x[0], 0) and on_boundary

# Outlet Surface
class Outlet(fenics.SubDomain):
	def inside(self, x, on_boundary):
		return fenics.near(x[0], L) and on_boundary

# Top Surface
def top(x, on_boundary):
	return fenics.near(x[1], Hf) and on_boundary

# Bottom Front Surface
def bottom_f(x, on_boundary):
	return fenics.near(x[1], 0) and on_boundary	

# Bottom Rear Surface
def bottom_r(x, on_boundary):
	return fenics.near(x[1], Hf-Hr) and on_boundary

# Ramp Surface
def ramp(x, on_boundary):
	return on_boundary and x[0]>=Lf and x[0] <= (L-Lr) and not fenics.near(x[1], Hf)
	

boundaries = fenics.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)

Inlet().mark(boundaries, 1)
Outlet().mark(boundaries, 2)

ds = fenics.Measure('ds', domain=mesh, subdomain_data=boundaries)

n = fenics.FacetNormal(mesh)

##########################################################################
bc_noslip_1 = fenics.DirichletBC(mixed_space.sub(0), fenics.Constant((0.0, 0.0)), top)
bc_noslip_2 = fenics.DirichletBC(mixed_space.sub(0), fenics.Constant((0.0, 0.0)), bottom_f)
bc_noslip_3 = fenics.DirichletBC(mixed_space.sub(0), fenics.Constant((0.0, 0.0)), ramp)
bc_noslip_4 = fenics.DirichletBC(mixed_space.sub(0), fenics.Constant((0.0, 0.0)), bottom_r)

bcs = [bc_noslip_1, bc_noslip_2, bc_noslip_3, bc_noslip_4]
##########################################################################
# Define Strain-Rate Tensor
def epsilon(u):
	grad_u = fenics.nabla_grad(u)
	return 0.5*(grad_u + grad_u.T)

##########################################################################
# Define Stress Tensor
def sigma(u, p):
	return 2*mu*epsilon(u) - p*fenics.Identity(len(u))

##########################################################################
# Construct Weak Forms
# Hydrodynamics
lhs = rho*fenics.dot(fenics.dot(un, fenics.nabla_grad(u)), v)*fenics.dx \
    + fenics.inner(2*mu*epsilon(u), epsilon(v))*fenics.dx \
	- p*fenics.div(v)*fenics.dx \
	+ fenics.div(u)*q*fenics.dx

rhs = fenics.dot(f,v)*fenics.dx \
	- p_left*fenics.dot(n,v)*ds(1) \
	- p_right*fenics.dot(n,v)*ds(2)

##########################################################################
iter = 0
iter_max = 10
converged = False
tol = 1E-8
eps = -1

# Solve Hydrodynamics w/ Picard Iteration
while not converged and iter < iter_max:
	iter += 1
	print('Iteration', iter)

	A = fenics.assemble(lhs)
	b = fenics.assemble(rhs)
	
	[bc.apply(A, b) for bc in bcs]

	fenics.solve(A, W.vector(), b)
	

	diff = W.sub(0).vector().get_local() - Wn.sub(0).vector().get_local()
	eps = numpy.linalg.norm(diff, ord=numpy.Inf)
	
	converged = eps < tol
	
	Wn.assign(W)
	print(diff)
1 Like
2

I believe the issue here is the use of Wn.split() instead of split(Wn). See the following simple example:

from dolfin import *

# Create function `u` in mixed space:
mesh = UnitSquareMesh(10,10)
VE = FiniteElement('CG', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,MixedElement(2*[VE,]))
u = Function(W)

# Two different ways to `split` the function:
u1, u2 = split(u)
u1s, u2s = u.split()

# Modify the original mixed function `u`:
u.assign(project(Constant((1,1)),W))

# Using `split(u)` reflects changes in `u`:
print(assemble(u1*dx))

# Using `u.split()` does not reflect changes:
print(assemble(u1s*dx))

This has been a known source of confusion for a while:

P.S. You can format code nicely by surrounding it with triple backticks and (optionally) specifying a language, as follows:

```python
line1
line2
etc.
```

Look up “Markdown” for more formatting tricks.

5 Likes
3

That did the trick! I’ve seen other places that the two splits did different things, but this helped clear it up. I’ve got it iterating, and even got a Newton iteration implemented. Now I’ve just got to get some kind of Anderson acceleration or other stabilization working! :stuck_out_tongue: