gistfile1.py

from dolfin import *
from dolfin_adjoint import *

# Create mesh
c = Point(0.0, 0.0, 0.0)
S = Sphere(c, 0.5)
mesh = Mesh(S, 10)

# Define mixed function spaces with Nedelec (edge) elements
X_h = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
X_N = X_h*X_h

# Set the control value
class Current(Expression):
    def eval(self, values, x):
        if abs(x[0])>DOLFIN_EPS and abs(x[1])>DOLFIN_EPS and abs(x[2])>DOLFIN_EPS:
            values[0] = -x[1]
            values[1] = x[0]
            values[2] = 0.0
        else:
            values[0] =0.0
            values[1] =0.0
            values[2] =0.0
    def value_shape(self):
        return(3,)

# solve the PDE
def solve_pde(fr, fi):
    sol = Function(X_N, name="Solution")
    (phi_hr, phi_hi) = TestFunctions(X_N)
    (E_hr, E_hi) = TrialFunctions(X_N)
    # Define variational form
    # real part
    F_r =  inner(curl(E_hr), curl(phi_hr)) * dx \
           - inner(curl(E_hi), curl(phi_hi)) * dx \
           + inner(E_hr, phi_hi) * dx \
           + inner(E_hi, phi_hr) * dx
    # imaginary part
    F_i = inner(curl(E_hr), curl(phi_hi)) * dx \
           + inner(curl(E_hi), curl(phi_hr)) * dx \
           - inner(E_hr, phi_hr) * dx \
           + inner(E_hi, phi_hi) * dx
    # left side
    a = F_r + F_i
    # right side
    L =  inner(fr, phi_hi) * dx + inner(fi, phi_hr) * dx - inner(fr, phi_hr) * dx + inner(fi, phi_hi) * dx
    solve(a == L, sol)
    return sol

# Optimization
def solve_optimal_control(n):
    # Instantiate the control value
    fr = Function(X_h)
    fi = Function(X_h)
    fr.assign(interpolate(Current(),X_h))
    fi.assign(interpolate(Current(),X_h))
    # Instantiate the left side
    y_re = Function(X_h)
    y_re.assign(Expression(('-x[1]', 'x[0]', '0.0')))
    y_ie = Function(X_h)
    y_ie.assign(Expression(('-x[1]', 'x[0]', '0.0')))
    #### Run the forward once to annotate
    sol = solve_pde(fr, fi)
    y_r, y_i = sol.split(deepcopy=True)
    #### The forward end
    Jr = Functional(0.5*inner(y_r-y_re, y_r-y_re)*dx + n/2 * inner(fr, fr) * dx)
    Ji = Functional(0.5*inner(y_i-y_ie, y_i-y_ie)*dx + n/2 * inner(fi, fi) * dx)
    m_r = TimeConstantParameter(fr)
    m_i = TimeConstantParameter(fi)
    Jm_r = assemble(0.5*inner(y_r-y_re, y_r-y_re)*dx + n/2 * inner(fr, fr) * dx)
    Jm_i = assemble(0.5*inner(y_i-y_ie, y_i-y_ie)*dx + n/2 * inner(fi, fi) * dx)
    dJdm_r = compute_gradient(Jr, m_r, forget=False)
    dJdm_i = compute_gradient(Ji, m_i, forget=False)
    #### Taylor test
    #def Jhat(m_r,m_i):
    #u = solve_pde(fr=m_r, fi=m_i)
    #u_r, u_i = u.split(deepcopy=True)
    #temp = 0.5*inner(u_r-y_re, u_r-y_re)*dx + n/2 * inner(fr, fr) * dx + 0.5*inner(u_i-y_ie, u_i-y_ie)*dx + n/2 * inner(fi, fi) * dx
    #return assemble(temp)

    #minconv = taylor_test(Jhat, (m_r,m_i), Jm_r, dJdm_r)
    #### Taylor test
    # Run the optimization
    reduced_functional = ReducedFunctional(Jr+Ji, [SteadyParameter(fr),SteadyParameter(fi)])
    f_opt = minimize(reduced_functional, method = 'Newton-CG', tol=1e-20, options={'maxiter': 200, 'disp': True})

    fr_opt, fi_opt = f_opt

    y_opt = solve_pde(fr=fr_opt, fi=fi_opt)


solve_optimal_control(1)