Tracking a curve of discontinuity of a function

Version:-
DOLFIN version: 2019.1.0 run using Conda

Motivation:-
I have a code where I have to evolve a discontinuous variable (\theta). At each step, I will evolve this variable via gradient descent,
\theta_{t+1}=\theta_t-\lambda\dfrac{\partial E^t}{\partial\theta},
where \lambda\approx0.01 is the learning rate.
The variable \theta is multivalued, but for my purposes, I want to restrict it to [0,2\pi). Note that the actual \theta being multivalued-continuous, the gradient is well defined. But restricting it to [0,2\pi), means that we get a “fake” gradient at the ‘line’ of discontinuity. Thus, I want to modify \dfrac{\partial E^t}{\partial\theta} near the discontinuity. This is the reason I need to keep track of this discontinuity.

Approach:-
My current approach was to read in the function, return the values at the nodes, and test which nodal values are 0 or 2\pi. However, this approach fails since it only returns 1 node.

location of 0's (array([840]),)
location of 2*pi's (array([], dtype=int64),)

Please check the MWE below.

import dolfin
print(f"DOLFIN version: {dolfin.__version__}")
from dolfin import *
import numpy as np


#Parameters
lx = float(1.0)
ly = float(1.0)
Nx = int(40)
Ny = int(40)
#c_r = float(0.1)
#Ref_No = int(input("Refinement number? -->"))


#Create mesh and define function space
mesh = RectangleMesh(Point(0., 0.), Point(lx, ly), Nx, Ny) 
x = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "Lagrange", 1)

##Reading in function
T = interpolate( Expression('atan2(-0.5*(x[0]-0.5*lx)+rt*0.5*(x[1]-0.5*ly), 0.5*(x[1]-0.5*ly)+rt*0.5*(x[0]-0.5*lx) )', rt=np.sqrt(3), lx=lx, ly=ly, degree=1), V)

###---------------------------------------------------------------------------------------------------------------
Tarray = T.vector()[:]


print("location of 0's", np.where( np.array(Tarray) == 0))
print("location of 2*pi's", np.where( np.array(Tarray) == 2*np.pi))                                                                     

Help Request:-
Any idea how i can determine the discontinuities? I’m stuck.

Why are you checking where it is exactly zero? Don’t you want to check where the function is outside the range? i.e. less than zero.
Consider the following:

import dolfin

print(f"DOLFIN version: {dolfin.__version__}")
from dolfin import *
import numpy as np


# Parameters
lx = float(1.0)
ly = float(1.0)
Nx = int(40)
Ny = int(40)
# c_r = float(0.1)
# Ref_No = int(input("Refinement number? -->"))


# Create mesh and define function space
mesh = RectangleMesh(Point(0.0, 0.0), Point(lx, ly), Nx, Ny)
x = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "Lagrange", 1)

##Reading in function
T = interpolate(
    Expression(
        "atan2(-0.5*(x[0]-0.5*lx)+rt*0.5*(x[1]-0.5*ly), 0.5*(x[1]-0.5*ly)+rt*0.5*(x[0]-0.5*lx) )",
        rt=np.sqrt(3),
        lx=lx,
        ly=ly,
        degree=1,
    ),
    V,
)

###---------------------------------------------------------------------------------------------------------------
Tarray = T.vector()[:]

try:
    import ufl
except ModuleNotFoundError:
    import ufl_legacy as ufl

uh = project(ufl.conditional(ufl.lt(T, 0), 1, 0), V)
print(uh.vector().get_local())


print("location of 0's", np.where(np.array(Tarray) == 0))
print("location of 2*pi's", np.where(np.array(Tarray) == 2 * np.pi))

File("T.pvd") << T
File("uh.pvd") << uh

Yielding

1 Like

Thank you. I was headed in the wrong direction.