(Deleted) Coding a primal dual active set algorithm


115
views
1
4 months ago by
b.g.  
Let \( \Omega = [0, 1]^2 \) and denote by \( \Omega_h \) the associated mesh. Let \( A^{-1}_p := \Omega_h \), \( A^{-1}_c \subset \Omega_h \), \(\psi := -0.075 \), \( c := 10^{-8} \), \( \gamma:= 0.011\) and \(\delta:= 0.01 \).

I am struggling in coding the following (primal-dual active set) algorithm:

Find \((u^n, p^n, \lambda^n)\) for \(n = 0, 1, \ldots \) such that:
\begin{align}
&R u^n - Mf + \pi(p^n) - \lambda^n = 0 \\
&p^n = \frac{\gamma}{\delta} \chi_{A_p^{n-1}} \\
&u^n = \psi \ \text{ on } \ A^{n-1}_c, \quad \lambda^n = 0 \ \text{ on } \ \Omega \setminus A^{n-1}_c \\
&\text{Update the sets:}\\
&A^n_c = \{x \in \Omega_h: \lambda^n(x) > c(u^n(x)-\psi) \} \\
&A^n_p = \{x \in \Omega_h: u^n(x) \leq \psi + \delta \}.
\end{align}

Here \( \chi_A(x) = 1 \) if  \( x \in A \) and  \( 0 \) otherwise,  \( R \) denotes the stiffness matrix for the (minus) Laplace operator,  \( M \) the mass matrix, and  \( (\pi(p^n))_i:= \int_\Omega \chi_{A^{n-1}_p}(x) \phi_i(x) dx \) , where  \( \{\phi_i\}_{i=1}^N \) is the finite element basis. I am using  \( P_1 \) Lagrange elements and homogeneous Dirichlet boundary conditions.

The strategy I implemented is the following: I loop over all of the mesh points and solve the first three equations pointwise. First, I compute  \( p^n \), then I fill  \( u^n \) and \( \lambda^n \) with the known values depending on where the mesh point lies, then finally I assemble the new matrix and load vector in order to solve the linear system in the first equation for the remaining unknowns.

At the moment I encounter the following problems.

1). My code yields a solution \( u < 0 \) on the boundary.

2). Unlike the paper from which I took this algorithm, the sets \( A_c \) and \( A_p \) stop changing values from the second or third iteration (it should normally take at least 6-7 iterations, and they should decrease).

I've only been using FEniCS for a week, so any hints or corrections are welcome. Thank you.

The code: 
from fenics import *
import numpy as np
from copy import deepcopy

gamma = 11*pow(10, -3); delta = pow(10, -2); psi = -0.075; c = pow(10, -8)

nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

n = V.dim()				
d = mesh.geometry().dim()		
dof_coord = V.tabulate_dof_coordinates()
dof_coord.resize((n, d))		

bc = DirichletBC(V, Constant(0.), lambda x, on_boundary: on_boundary)

# Variational problem
u = TrialFunction(V)
v = TestFunction(V)
p = Function(V)				
lamda = Function(V)		

a = dot(grad(u), grad(v))*dx
R = assemble(a)
M = assemble(u*v*dx)
bc.apply(R)
bc.apply(M)
M = M.array()				# Mass matrix
R = R.array()				# Stifness matrix

Ap = deepcopy(dof_coord)
Ac = np.array([x if x[0] in {.375, .5, .625} and x[1] in {.375, .5, .625} else None for i, x in enumerate(dof_coord)])

# Having taken A_p^{-1} := Omega, we may set
p.vector()[:] = gamma/delta*np.ones(n)

u = Function(V)				
vtkfile = File('pdas.pvd')		

index_u = []; index_lamda = []		
R_ = np.zeros((n,n)); b_ = np.zeros(n)  

max_iter = 4
for iter in range(max_iter):
	Pi = assemble(p*v*dx); bc.apply(Pi); Pi = Pi.get_local()
	index_lamda = [i for i, x in enumerate(dof_coord) if any((x==y).all() for y in Ac)]
	index_u = list(set(range(n))-set(index_lamda))
	for i, x in enumerate(dof_coord):
		p.vector()[i] = gamma/delta if any((x==y).all() for y in Ap) else 0.
		if any((x==y).all() for y in Ac):
			u.vector()[i] = psi
		else:
			lamda.vector()[i] = 0.0
		b_[i] = -1*M[i].dot(np.ones(n))-Pi[i]-psi*np.sum([R[i][j] for j in index_lamda])
		for j in range(n):
			if j in index_u:
				R_[i][j] = R[i][j]
			elif j in index_lamda and j == i:
				R_[i][j] = -1
			else:
				R_[i][j] = 0
	# Compute the remaining unknowns for u and lambda
	U = np.linalg.solve(R_, b_)
	for l, z in enumerate(U):
		if l in index_u:
			u.vector()[l] = z
		else:
			lamda.vector()[l] = z
	index_u, index_lamda = [], []
	# Update the active sets
	Ap = np.array([x if u(x) <= psi+delta else None for x in dof_coord])
	Ac = np.array([x if lamda(x) > c*(u(x)-psi) else None for x in dof_coord])
	vtkfile << u​

Community: FEniCS Project
Please login to add an answer/comment or follow this question.
The thread is closed. No new answer/comment may be added.

Similar posts:
Search »
  • Nothing matches yet.