Commit e8936c0f authored by Élise Le Mélédo's avatar Élise Le Mélédo
Browse files

New version: Added redistancing, corrected major bug for very high order,...

New version: Added redistancing, corrected major bug for very high order, updated documentation, added test cases. Warning: use the new meshes from this version on.
parent 521259b8
Pipeline #4230 passed with stage
in 10 seconds
......@@ -5,7 +5,7 @@
# ==============================================================================
# Pythonic interfacing constraints (safety barrier)
# ==============================================================================
__all__ = ["BarycentricBasis","BarycentricGradients"]
__all__ = ["BarycentricBasis", "BarycentricGradients"]
# ==============================================================================
# Preliminary imports
......@@ -63,9 +63,6 @@ def BarycentricGradients(Nele, Type, Order, points, n, vol):
else: raise(ValueError("Error: Barycentric basis functions are not yet implemented for elements having "+str(Nele)+"edges. Abortion."))
################################################################################
################################################################################
# ==============================================================================
# Technical routines for the basis functions computations
# ==============================================================================
......@@ -93,6 +90,7 @@ def SimplicialBasis(Type, Order, points):
# Generate the basis functions according to the wished order (automatically)
if Type == "B":
# ------ Generate the permutations
# Generate the indices of the vertices involved in the generation of the barycentric coordinates
BaryIndex = np.arange(3)
# Generate the combinations corresponding to each barycentric coordinate to select in the monomial generation
......@@ -100,10 +98,34 @@ def SimplicialBasis(Type, Order, points):
# Initilalising the basis function according to the number of basis functions
EvaluatedBase = np.zeros((points.shape[0], len(ConsideredIndices)))
# ------ Compute the basis functions
# Compute the basis function in their general form
for (i,Cind) in enumerate(ConsideredIndices):
EvaluatedBase[:, i] = (np.math.factorial(Order)/np.prod([np.math.factorial(Cind.count(e)) for e in set(Cind)]))*np.prod(points[:,Cind], axis=1)
# ------ Reorder the basis functions
# Map the functions to match the Barycentric coordinates initiated by the MeshStructure (Laborious but working)
# Find the order on each edge
Reordering = np.zeros(len(ConsideredIndices), dtype=int); Reordering[0] = 1
Tmp1 = np.zeros(Order+1); Tmp1[0] = 1;
for j in range(1, Order+1): Tmp1[j] = Tmp1[j-1] + j
Tmp2 = np.zeros(Order); Tmp2[0] = 1;
for j in range(1, Order): Tmp2[j] = Tmp2[j-1] + j+1;
# Assigning the reordering associated to the Meshstrcture
Reordering[0:Order+1] = Tmp1
Reordering[Order+1:2*Order+1] = np.arange(Reordering[Order]+1, len(ConsideredIndices)+1)
Reordering[2*Order+1:3*Order] = Tmp2[::-1][:-1]
# Assigning the reordering associated to the internal basis functions
Index = np.setdiff1d(np.arange(len(ConsideredIndices))+1,Reordering)
Reordering[int(3*Order):] = Index
Reordering = Reordering-1
# Reorder the basis function accordingly
EvaluatedBase = EvaluatedBase[:, Reordering]
# Safety check
else: raise NotImplementedError("The type of wished basis function is not implemented yet. Abortion.")
......@@ -144,6 +166,7 @@ def SimplicialGradients(Type, Order, points, n, vol):
# Generate the basis functions according to the wished order (automatically)
if Type == "B":
# ------ Generate the permutations
# Generate the indices of the vertices involved in the generation of the barycentric coordinates
BaryIndex = np.arange(3)
# Generate the combinations corresponding to each barycentric coordinate to select in the monomial generation
......@@ -152,6 +175,7 @@ def SimplicialGradients(Type, Order, points, n, vol):
EvaluatedGradient = np.zeros((points.shape[0], len(ConsideredIndices), 3))
Grad = np.zeros((points.shape[0], 2, len(ConsideredIndices)))
# ------ Compute the components involved in each gradient
# Compute the derivatives of the basis function in their general form
for (j,Cind) in enumerate(ConsideredIndices):
for k in range(3):
......@@ -160,6 +184,30 @@ def SimplicialGradients(Type, Order, points, n, vol):
Ccind.remove(k)
EvaluatedGradient[:, j, k] = np.array([(np.math.factorial(Order)/np.prod([np.math.factorial(Ccind.count(e)) for e in set(Ccind)]))*np.prod(points[:,Ccind], axis=1)])
# ------ Reorder the gradients functions
# Map the functions to match the Barycentric coordinates initiated by the MeshStructure (Laborious but working)
# Find the order on each edge
Reordering = np.zeros(len(ConsideredIndices), dtype=int); Reordering[0] = 1
Tmp1 = np.zeros(Order+1); Tmp1[0] = 1;
for j in range(1, Order+1): Tmp1[j] = Tmp1[j-1] + j
Tmp2 = np.zeros(Order); Tmp2[0] = 1;
for j in range(1, Order): Tmp2[j] = Tmp2[j-1] + j+1;
# Assigning the reordering associated to the Meshstrcture
Reordering[0:Order+1] = Tmp1
Reordering[Order+1:2*Order+1] = np.arange(Reordering[Order]+1, len(ConsideredIndices)+1)
Reordering[2*Order+1:3*Order] = Tmp2[::-1][:-1]
# Assigning the reordering associated to the internal basis functions
Index = np.setdiff1d(np.arange(len(ConsideredIndices))+1,Reordering)
Reordering[int(3*Order):] = Index
Reordering = Reordering-1
# Reorder the basis function accordingly
EvaluatedGradient = EvaluatedGradient[:, Reordering,:]
# ------ Adapt the gradients to the geometry
# Scale the values to the geometry
Grad = np.dot(np.swapaxes(np.array([EvaluatedGradient[:,:,0]-EvaluatedGradient[:,:,2],
EvaluatedGradient[:,:,1]-EvaluatedGradient[:,:,2]]),0,2),\
......
"""
Implementation of the CoCaIn-BPG algorithm (Convex-Concave Backtracking for Inertial Bregman Proximal Gradient Algorithms),
suitable for non-convex problem, without proximal mapping.
Code adapted from the example given by Mahesh Chandra Mukkamala, Peter Ochs, Thomas Pock and Shoham Sabach and related
to their paper that can be found on https://arxiv.org/abs/1904.03537
.. note::
The global optimiser may not work if initial guess is saddle point of local extrema
The algoritm will converge to a critical point (where Grad(f)=0)
"""
# ==================================================================================
# Preliminary imports
# ==================================================================================
# Classical library imports
import numpy as np
# ==================================================================================
# Solver definition
# ==================================================================================
class Optimiser():
def __init__(self):
self.del_val = 0.15
self.eps_val = 0.00001
self.uL_est = 1/(1-self.del_val) + 1e-3
self.lL_est = 1e-4*self.uL_est
self.max_iter = 100
# -----------------------------------------------------------------
# Helper functions
# -----------------------------------------------------------------
def breg(self, U, U1):
temp = 0.5*(np.sum(np.multiply(U-U1,U-U1)))
if temp >=1e-15: return temp
else: return 0
def abs_fun(self, fun, grad, U, U1):
G = grad(U1)
x,y = U[0], U[1]
x_1, y_1 = U1[0], U1[1]
return fun(U1)+np.sum(np.multiply(G,U-U1))
def make_update(self, y, evalgrad, uL_est):
temp_p = y - (1/uL_est)*evalgrad
return temp_p
def find_gamma(self, U, prev_U, uL_est, lL_est):
gamma = 1
kappa = (self.del_val - self.eps_val)*(uL_est/(uL_est+lL_est))
y_U = U + gamma*(U-prev_U) # extrapolation
while (kappa*self.breg(prev_U, U)< self.breg(U, y_U)):
gamma = gamma*0.9
y_U = U + gamma*(U-prev_U)
return y_U, gamma
def do_lb_search(self, fun, grad, U, U1, uL_est, lL_est):
y_U, gamma = self.find_gamma(U, U1, uL_est, lL_est)
while((self.abs_fun(fun, grad, U, y_U) - fun(U) - (lL_est*self.breg(U, y_U))) > 1e-7):
lL_est = 2*lL_est
y_U, gamma = self.find_gamma(U, U1, uL_est, lL_est)
return lL_est, y_U, gamma
def do_ub_search(self, fun, grad, y_U, uL_est):
grad_u = np.squeeze(grad(y_U))
x_U = self.make_update(y_U, grad_u, uL_est)
while((self.abs_fun(fun, grad, x_U, y_U)- fun(x_U)+ (uL_est*self.breg(x_U, y_U))) < -1e-7):
uL_est = (2)*uL_est
x_U = self.make_update(y_U, grad_u, uL_est)
return uL_est, x_U
# -----------------------------------------------------------------
# Main optimisation routine
# -----------------------------------------------------------------
def Optimise(self, fun, grad, U0, Point, radius):
uL_est = self.uL_est
lL_est = self.lL_est
U = U0
init_U = U
prev_U = U
tmpU = U+1
gamma_vals = [0]
uL_est_vals = [uL_est]
lL_est_vals = [uL_est]
temp = fun(U)
func_vals = [temp]
lyapunov_vals = [temp]
U_vals = [init_U]
k=0
while np.any(abs(U-tmpU)>1e-8) and k<10:
k=k+1
lL_est, y_U, gamma = self.do_lb_search(fun, grad, U, prev_U, uL_est, lL_est)
prev_U = U
tmpU=U
uL_est, U = self.do_ub_search(fun, grad, y_U, uL_est)
temp = fun(U)
uL_est_vals = uL_est_vals + [uL_est]
lL_est_vals = lL_est_vals + [lL_est]
gamma_vals = gamma_vals + [gamma]
U_vals = U_vals + [U]
if np.isnan(temp): return(prev_U)
InnerCircle = np.linalg.norm(U-Point)
if InnerCircle<=radius: U= U
else: U = Point+radius*(U-Point)/InnerCircle
func_vals = func_vals + [temp]
lyapunov_vals = lyapunov_vals + [(1/uL_est)*temp + self.del_val*self.breg(U, prev_U)]
return(U)
"""
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
|
Module that computes the solution of any equation of the type \frac{\partial \phi}{\partial t} + H(\phi, t) = 0 over
a given Cartesian domain.
It implements the WENO scheme for Cartesian grids for HJB equations, following the work of Peng and Jian,
"Weighted ENO schemes for HJB equations, PII: S106482759732455X", and uses a simple forward euler method in time.
"""
# ==============================================================================
# Preliminary imports
# ==============================================================================
import copy
import numpy as np
import matplotlib.pyplot as plt
# ==============================================================================
# WENO tools
# ==============================================================================
# --- Interate the equation (7) until the BregmanInitTime by the Weno 5th order scheme for the residuals, and Euler explicit in time (one iteration)
# Defining a consistant flux with H, here Global LFX
def Hh(dxP, dxM, dyP, dyM, H, Hx, Hy):
alph = np.max([np.max(Hx(dxP, dyP)), np.max(Hx(dxP, dyM)), np.max(Hx(dxM, dyP)), np.max(Hx(dxM, dyM))])
bet = np.max([np.max(Hx(dxP, dyP)), np.max(Hx(dxP, dyM)), np.max(Hx(dxM, dyP)), np.max(Hx(dxM, dyM))])
Val = H(0.5*(dxP+dxM), 0.5*(dyP+dyM)) - alph*(0.5*(dxP-dxM)) - bet*(0.5*(dyP-dyM))
return(Val)
# Defining the Weno approximation function
def W(a,b,c,d, eps):
alph0 = 1/((eps+13*(a-b)**2+3*(a-3*b)**2)**2)
alph1 = 6/((eps+13*(b-c)**2+3*(b+c)**2)**2)
alph2 = 3/((eps+13*(c-d)**2+3*(3*c-d)**2)**2)
w0 = alph0/(alph0+alph1+alph2)
w2 = alph2/(alph0+alph1+alph2)
return(1/3*w0*(a-2*b+c)+1/6*(w2-0.5)*(b-2*c+d))
# ==============================================================================
# Actual routine
# ==============================================================================
def Iterate(xx, yy, InitVal, Tmax, cfl, eps, H, Hx, Hy):
# Retrieving the grid size
DeltaX = np.abs(xx[1]-xx[0])
DeltaY = np.abs(yy[1]-yy[0])
# Initialising the iterated value
Val = copy.deepcopy(InitVal[:,:])
t=0
while t<Tmax:
# Registering the buffer term to impose boundary conditions
Buff = copy.deepcopy(Val[:,:])
# Computing the differences
DeltaxPPhi = np.zeros(Val.shape)
DeltaxMPhi = np.zeros(Val.shape)
DeltayPPhi = np.zeros(Val.shape)
DeltayMPhi = np.zeros(Val.shape)
DeltaxPPhi[:,:-1] = Val[:, 1:] - Val[:, :-1]
DeltaxMPhi[:, 1:] = Val[:, 1:] - Val[:, :-1]
DeltayPPhi[:-1,:] = Val[1:,:] - Val[:-1, :]
DeltayMPhi[1:, :] = Val[1:,:] - Val[:-1, :]
DDeltaxMPPhi = np.zeros(InitVal.shape)
DDeltayMPPhi = np.zeros(InitVal.shape)
DDeltaxMPPhi[:,1:] = DeltaxPPhi[:,1:]-DeltaxPPhi[:,:-1]
DDeltayMPPhi[1:,:] = DeltayPPhi[1:,:]-DeltayPPhi[:-1,:]
# Approximating the derivatives
dxPphi = np.zeros(InitVal.shape)
dyPphi = np.zeros(InitVal.shape)
dxMphi = np.zeros(InitVal.shape)
dyMphi = np.zeros(InitVal.shape)
dxPphi[:,2:-2] = 1/(12*DeltaX)*(-DeltaxPPhi[:,:-4]+7*DeltaxPPhi[:,1:-3]+7*DeltaxPPhi[:,2:-2]-DeltaxPPhi[:,3:-1])\
+W(DDeltaxMPPhi[:,4:]/DeltaX, DDeltaxMPPhi[:,3:-1]/DeltaX, DDeltaxMPPhi[:,2:-2]/DeltaX, DDeltaxMPPhi[:,1:-3]/DeltaX, eps)
dyPphi[2:-2,:] = 1/(12*DeltaY)*(-DeltayPPhi[:-4,:]+7*DeltayPPhi[1:-3,:]+7*DeltayPPhi[2:-2,:]-DeltayPPhi[3:-1,:])\
+W(DDeltayMPPhi[4:,:]/DeltaY, DDeltayMPPhi[3:-1,:]/DeltaY, DDeltayMPPhi[2:-2,:]/DeltaY, DDeltayMPPhi[1:-3,:]/DeltaY, eps)
dxMphi[:,2:-2] = 1/(12*DeltaX)*(-DeltaxPPhi[:,:-4]+7*DeltaxPPhi[:,1:-3]+7*DeltaxPPhi[:,2:-2]-DeltaxPPhi[:,3:-1])\
-W(DDeltaxMPPhi[:,:-4]/DeltaX, DDeltaxMPPhi[:,1:-3]/DeltaX, DDeltaxMPPhi[:,2:-2]/DeltaX, DDeltaxMPPhi[:,3:-1]/DeltaX, eps)
dyMphi[2:-2,:] = 1/(12*DeltaY)*(-DeltayPPhi[:-4,:]+7*DeltayPPhi[1:-3,:]+7*DeltayPPhi[2:-2,:]-DeltayPPhi[3:-1,:])\
-W(DDeltayMPPhi[:-4,:]/DeltaY, DDeltayMPPhi[1:-3,:]/DeltaY, DDeltayMPPhi[2:-2,:]/DeltaY, DDeltayMPPhi[3:-1,:]/DeltaY, eps)
# Compute the time step due to cfl
lbd = np.max([np.max(Hx(dxPphi[2:-2,:] , dyPphi[2:-2,:] )), np.max(Hx(dxPphi[2:-2,:] , dyMphi[2:-2,:] )), np.max(Hx(dxMphi[2:-2,:] , dyPphi[2:-2,:] )), np.max(Hx(dxMphi[2:-2,:] , dyMphi[2:-2,:] ))])
dt = cfl*np.max([DeltaX,DeltaY])/lbd
if (Tmax<t+dt): dt = Tmax-t
# Computing the residuals and updating with straighforward Euler
dphi = -Hh(dxPphi[2:-2, 2:-2], dxMphi[2:-2, 2:-2], dyPphi[2:-2, 2:-2], dyMphi[2:-2, 2:-2], H, Hx, Hy)
Val[2:-2,2:-2] = Val[2:-2,2:-2] + dt*dphi
Val[0:2,:] = Buff[0:2,:]
Val[-2:,:] = Buff[-2:,:]
Val[:,0:2] = Buff[:,0:2]
Val[:,-2:] = Buff[:,-2:]
t = t+dt
# Retrieve the gradients of the approximation at BregmanInitTime at the grid points
Val[0:2,:] = 0
Val[-2:,:] =0
Val[:,0:2] = 0
Val[:,-2:] = 0
dxVal = 0.5*(dxPphi+dxMphi)
dyVal = 0.5*(dyPphi+dyMphi)
# Plot the redistanced level set
return(Val, dxVal, dyVal)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment