extract dual potentials (φ, ψ) and their derivatives using POT #728

I've been using POT to solve 2D optimal transport problems involving transforming shapes (ellipse to circle, etc.), and I’ve run into an issue I couldn’t resolve after exploring several of the available APIs.
I tried ot.sinkhorn and the convolutional barycenter. however for the first didn't return reliable values, and the second seems to only return the barycenter.

This is an example of how I could define my source and target:
def S1(x, y, r):
O = np.zeros(len(x))
for k in range(len(x)):
if x[k]**2 + y[k]2 < r2:
O[k] = 1
return O / np.sum(O)

def S2(x, y, a, b):
O = np.zeros(len(x))
for k in range(len(x)):
if (x[k]2) / a2 + (y[k]2) / b2 < 1:
O[k] = 1
return O / np.sum(O)

r = 0.5
a = 0.5
b = 0.3

x = np.linspace(-1, 1, 50)
X_p, Y_p = np.meshgrid(x, x)
X_p_flat = X_p.flatten()
Y_p_flat = Y_p.flatten()

mu_s = S1(X_p_flat, Y_p_flat, r)
mu_t = S2(X_p_flat, Y_p_flat, a, b)

What I'm struggling with is a clear, stable example showing how to extract φ and ∇φ from any of the above APIs, or others.