### (Deleted) Coding a primal dual active set algorithm

115
views
1
4 months ago by
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.