Using Dolfin-adjoint for optimization problem without PDE constraint

dolfin-adjoint
1

Hello all,
I am trying to use Dolfin-adjoint to run an optimization problem that does not have any PDE constraint. Below is my code:

import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from dolfin_adjoint import *
import moola

n = 10
mesh = RectangleMesh(Point(-1,-1),Point(1,1), n, n)

V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
S0 = Constant(1)

bc = DirichletBC(V, 1, "on_boundary")

J = assemble((0.5*inner(grad(u), grad(u)) - u*S0)*dx)
J_hat = ReducedFunctional(J, Control(u))   
u_opt = minimize(J_hat, method = "L-BFGS-B", options = {"gtol": 1e-16, "ftol": 1e-16})

u_opt.rename("u_opt","temp")
fileY = File("temp.pvd");
fileY << u_opt;

However, I don’t know how I should apply the boundary condition of the control variable “u” to the minimization problem. Also, I am pretty sure that the above code is not giving me the correct answer because when I am trying to solve its Euler-Lagrange equation using the below code, the answers are different:

from dolfin import *
from fenics import *
import numpy as np

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

mesh = RectangleMesh(Point(-1,-1),Point(1,1), 10, 10)
V = FunctionSpace(mesh, "CG", 1)

du = TrialFunction(V)            # Incremental displacement
u  = Function(V) 
v  = TestFunction(V)             # Test function
S0 = Constant(1)

F =  (inner(grad(u), grad(v)) - v*S0)*dx
bc = DirichletBC(V, 1, "on_boundary")
solve(F == 0, u, bc)
 
u.rename("u","temp")
fileY = File("temp.pvd");
fileY << u;

Thanks in advance for your help and support!

1 Like
2

This might not be the most elegant approach, but I would suggest using a immediate variable between your control and the function you would like to assemble over, as shown in the following example:

import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from dolfin_adjoint import *
import moola

n = 10
mesh = RectangleMesh(Point(-1,-1),Point(1,1), n, n)

V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
S0 = Constant(1)

bc = DirichletBC(V, 1, "on_boundary")

v = project(u, V, bcs=bc)
J = assemble((0.5*inner(grad(v), grad(v)) - v*S0)*dx)
J_hat = ReducedFunctional(J, Control(u))   
u_opt = minimize(J_hat, method = "L-BFGS-B", options = {"gtol": 1e-16, "ftol": 1e-16})
J_hat(u_opt)
fileY = File("temp.pvd");
fileY << v.block_variable.saved_output;
print(assemble(v*ds))
1 Like
3

Thank you very much for your very fast and great help!