I have a few parameters in an Expression that I would like to use as control variables. For example,
from dolfin import *
from dolfin_adjoint import *
c = Constant(1)
coeff = Expression('c*cos(x[0])', c=c, degree=1)
m = Control(coeff.c)
but this doesn’t work. How can we implement something like this? Many thanks in advance.
Hi,
When you use an expression, You have to supply the derivative of that expression. (also, the control should be m=Control(c)) See for instance:
http://www.dolfin-adjoint.org/en/release/documentation/time-dependent-wave/time-dependent-wave.html
However, at least in your simple case, you can avoid using Expression by using SpatialCoordinate:
from dolfin import *
from dolfin_adjoint import *
mesh = UnitSquareMesh(10,10)
c = Constant(1)
x = SpatialCoordinate(mesh)
coeff = c*cos(x[0])
m = Control(c)
Jhat = ReducedFunctional(assemble(coeff**2*dx(domain=mesh)), m)
dJdm = Jhat.derivative()
Hi @dokken,
Thanks for the link to the example. My expressions are indeed more complicated than my MWE above, but link pointed me in the right direction. Specifically, my problem involves vector valued Expressions, as well as the spatial derivatives of those expressions to describe a parametric surface whose cross section is the arc of circle. I have been using sympy to generate the expressions and their derivatives.
On a related note, in the code below, I have an implicit relationship between variables W, L and \eta=, currently approximated by a fitted curve. Is there any trick out to wrap this implicit relationship into an Expression object? (Because this is a hack and hard to generalize.) As I mention in the code comments, I can use sympy nsolve to get the value of \eta for a given w and L running the model in “forward only mode”, but I don’t see how this would help me with defining the derivatives for use with dolfin adjoint.
Thanks again!
from dolfin import *
from dolfin_adjoint import *
import sympy as sp
import numpy as np
DEGREE = 3
def ccode(z):
return sp.printing.ccode(z)
class AdjointGeometry(object):
def __init__(self, w=None, L=None):
w_val = w # given value, does not change
self.length = L # this will be used as a control variable
W, eta, L, x1, x2 = sp.symbols('w, eta, length, x[0], x[1]')
CC = W/L # width to length ratio
# constant \eta is given by the following implicit relationship
# -CC * \eta = \sqrt(2*(1-cos(\eta)))
# previously, in "forward-only mode" I have been able to provide numerical value for use sympy
to solve as such:
# eta = sp.nsolve(-c*eta - 2*(1-sp.cos(eta)), eta, -5)
# For adjoint use, as temporary workaround, fitting explicit curve to desired portion as an
approximation:
from scipy import optimize
def test_func(x, a, b,c,d,e,f):
return a + b*x + c*x**2 + d*x**3 + e*np.cos(f*x)
a, b,c,d,e,f = sp.symbols('a, b,c,d,e,f')
eta = a + b*CC + c*CC**2 + d*CC**3 + e*sp.cos(f*CC)
eta_fit = np.linspace(-6.28,0.001,100)
x_data = np.sqrt(2*(1-np.cos(eta_fit))/eta_fit**2)
params, params_covariance = optimize.curve_fit(test_func, x_data, eta_fit, maxfev=10000)
# Now build the desired expressions (gamma and its spatial derivatives)
A = -0.5*W
B = 0.5*W * (sp.sin(eta)/(1-sp.cos(eta)))
R = L/eta
# store the symbolic components in a list for further use
gamma_sp = [None, None, None]
gamma_sp[0] = A*sp.cos(x1*L/R) + B*sp.sin(x1*L/R)
gamma_sp[1] = x2
gamma_sp[2] = A*sp.sin(x1*L/R) - B*sp.cos(x1*L/R)
self.gamma = Expression((ccode(gamma_sp[0]).replace('M_PI', 'pi'),
ccode(gamma_sp[1]).replace('M_PI', 'pi'),
ccode(gamma_sp[2]).replace('M_PI', 'pi')),
pi=pi,
length=self.length,
w=w_val,
a=params[0],
b=params[1],
c=params[2],
d=params[3],
e=params[4],
f=params[5],
degree=DEGREE,
name="gamma")
# spatial derivative of above gamma (needed for PDE), similarly for \partial\gamma / \partial x2
self.Gsub1 = Expression((ccode(gamma_sp[0].diff(x1)).replace('M_PI', 'pi'),
ccode(gamma_sp[1].diff(x1)).replace('M_PI', 'pi'),
ccode(gamma_sp[2].diff(x1)).replace('M_PI', 'pi')),
pi=pi,
length=self.length,
w=w_val,
a=params[0],
b=params[1],
c=params[2],
d=params[3],
e=params[4],
f=params[5],
degree=DEGREE,
name="Gsub1")
# Now need the derivatives of the above wrt length for dolfin-adjoint
self.d_gamma = Expression((ccode(gamma_sp[0].diff(L)).replace('M_PI', 'pi'),
ccode(gamma_sp[1].diff(L)).replace('M_PI', 'pi'),
ccode(gamma_sp[2].diff(L)).replace('M_PI', 'pi')),
pi=pi,
length=self.length,
w=w_val,
a=params[0],
b=params[1],
c=params[2],
d=params[3],
e=params[4],
f=params[5],
degree=DEGREE,
name="d_gamma")
self.d_Gsub1 = Expression((ccode(gamma_sp[0].diff(x1, L)).replace('M_PI', 'pi'),
ccode(gamma_sp[1].diff(x1, L)).replace('M_PI', 'pi'),
ccode(gamma_sp[2].diff(x1, L)).replace('M_PI', 'pi')),
pi=pi,
length=self.length,
w=w_val,
a=params[0],
b=params[1],
c=params[2],
d=params[3],
e=params[4],
f=params[5],
degree=DEGREE,
name="d_Gsub1")
# per the example, store the derivatives as dictionaries
self.gamma.dependencies = [self.length]
self.gamma.user_defined_derivatives = {self.length: self.d_gamma}
self.Gsub1.dependencies = [self.length]
self.Gsub1.user_defined_derivatives = {self.length: self.d_Gsub1}
# generate initial parametric geometry
w = Constant(0.5, name="width")
L = Constant(0.25*pi, name="length")
geo = AdjointGeometry(w=0.5, L=L)
m = Control(L)