Hi everyone,
I would like to compute the adjoint using dolfin-adjoint for forward solve which uses a higher order runge-kutta method. However, if I run the taylor_test using RK2, I get clearly incorrect convergence orders:
RK2 Taylor [3.7722606893233466, 0.3183533286949072, 1.5733204186895104].
For the Euler timestepping (=RK1) I seem to get the correct taylor test convergence orders:
RK1 Taylor [2.2290340416041228, 2.1297747367037365, 2.069543839265416]
What I checked so far:
- As far as I can tell, the RK2 is correct, I checked this by computing the convergence order which should be 2 for RK order 2 (test included in code below):
[1.99523983949585, 1.998757686994705, 1.9996831836595186, 1.9999203012013684, 1.999980163477117] - I also compared the computed gradient DJ for Euler (RK1) and RK2. The values seem to be the same to be the same. For DG order 2, m_size=[40,80]
- RK1: DJ =
[3.0960247841905786, 3.1187228838225205] - RK2: DJ =
[3.13882278636299, 3.140900383048848]
- RK1: DJ =
I really do not know what could cause the taylor test not to converge. I also created the dependency graphs with tape.visualize() and they seem to be correct for RK2. The control variable seems to be connected the right way to the weak formulation. (The graph is a bit too big to include here)
Any ideas on what could be wrong would be greatly appreciated ! 
I shorted my code as far as possible. I use the fenics+dolfin-adjoint 2019.01.0 version.
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
from dolfin_adjoint import Constant as Constant
import numpy as np
parameters["ghost_mode"] = "shared_facet"
parameters["form_compiler"]["cpp_optimize"] = True
class gD_derivative(UserExpression):
def __init__(self, gD, **kwargs):
""" Construct the source function derivative """
super().__init__(**kwargs)
self.t = 0.0
self.beta_x = gD.beta_x
self.gD = gD # needed to get the matching time instant
def eval(self, value, x):
""" Evaluate the source function's derivative """
t = self.gD.t
value[0] = -2 * pi * t * cos(2 * pi * (x[0] - self.beta_x * t))
def value_shape(self):
return ()
def forward(beta, gD, ic, f, unit_mesh, h, dg_order, time_stepping,
annotate=False, show_plots=False, show_final_plots=False):
dx = Measure('dx', domain=unit_mesh) # area
dS = Measure('dS', domain=unit_mesh) # interior boundaries
ds = Measure('ds', domain=unit_mesh) # exterior boundaries
V = FunctionSpace(unit_mesh, 'DG', dg_order)
u = TrialFunction(V)
v = TestFunction(V)
face_normal = FacetNormal(unit_mesh)
factor = 0.2
maxvel = beta[0]
print(float(beta.beta_x))
delta_t = (1 / (dg_order ** 2 + 1)) * h * factor * (1 / abs(float(beta.beta_x)))
def F_operator(beta, u_old, gD, v, f, n): # beta, u1, gD, v, f, face_normal
flux = beta * u_old # physical flux
F = dot(grad(v), flux) * dx + inner(v, f) * dx
flux_dS_ = dot(avg(flux), n('+')) + 0.5 * maxvel('+') * jump(u_old)
flux_ds_ = dot(0.5 * (flux + gD * beta), n) + 0.5 * maxvel * (u_old - gD)
F -= inner(flux_dS_, jump(v)) * dS
F -= inner(flux_ds_, v) * ds
return F
if time_stepping == 'RK2':
a = inner(v, u) * dx
u1 = project(gD, V, annotate=annotate)
u2 = Function(V, annotate=annotate)
L1 = F_operator(beta, u1, gD, v, f, face_normal) # u_old, beta, gD, v, f, n
L2 = F_operator(beta, u2, gD, v, f, face_normal)
else:
u1 = project(ic, V, annotate=annotate)
u1.rename('u1', 'u1')
const_delta_t = Constant(1 / delta_t, name='delta_t')
F = F_operator(beta, u, gD, v, f, face_normal) - inner(const_delta_t * (u - u1), v) * dx
a, L = lhs(F), rhs(F)
u_new = Function(V) # storage for computed time step
tend = 1.0 # set up time stepping
steps = int(tend / delta_t)
t = 0.0
gD.t = 0.0
for n in range(0, steps):
if n == steps - 1: # final step
delta_t = tend - t
if time_stepping == 'RK2': # runge-kutta 2
# pause_annotation()
solve(a == L1, u_new, annotate=annotate)
u2.assign(u1 + delta_t * u_new, annotate=annotate)
gD.t = gD.t + delta_t
# continue_annotation()
solve(a == L2, u_new, annotate=annotate)
u1.assign(0.5 * u1 + 0.5 * (u2 + delta_t * u_new), annotate=annotate)
else: # explicit euler
solve(a == L, u_new, annotate=annotate)
u1.assign(u_new, annotate=annotate)
t += delta_t
gD.t = t
if n % int(0.1 * steps) == 0:
print("step n=", n, "of steps=", steps)
if show_plots:
plot(project(gD, V), label='gD')
plot(u1, label='u1', linestyle='--')
plt.legend()
plt.show()
if show_final_plots:
plot(project(gD, V), label='gD')
plot(u1, label='u1', linestyle='--')
plt.legend()
plt.show()
print("t final", t, "delta_t", delta_t)
return u1, t
def set_up_example(m_size):
unit_mesh_ = UnitIntervalMesh(m_size)
source_term = Expression('0.0', degree=3)
ic = Expression('sin(2*pi*(x[0]))', degree=3, name="ic")
# constant flow speed
beta_x = Constant(0.25, name="beta_x")
beta = Expression(("beta_x",), degree=1, beta_x=beta_x, name="beta_vec") # as_vector((beta_x,))
beta.dependencies = [beta_x]
beta.user_defined_derivatives = {beta_x: Constant((1.0,))}
# setup boundary function
gD = Expression("sin(2 * pi * (x[0] - beta_x * t))", degree=1, t=0, beta_x=beta_x)
gD.dependencies = [beta_x]
gD.user_defined_derivatives = {beta_x: Expression("-2 * pi * t * cos(2 * pi * (x[0] - beta_x * t))",
degree=1, t=0, beta_x=beta_x)}
return unit_mesh_, source_term, ic, beta, beta_x, gD
def dolfin_adjoint_taylor_test():
dg_ord = 2
DJ = {'RK1': [], 'RK2': []}
rk_ord = 'RK1'
m_size = 20
unit_mesh_, source_term, ic, beta, beta_x, gD = set_up_example(m_size)
control = Control(beta_x)
u_final, t_final = forward(beta, gD, ic, source_term, unit_mesh_, h=1 / m_size, annotate=True,
dg_order=dg_ord, time_stepping=rk_ord, show_plots=False, show_final_plots=True)
error = errornorm(gD, u_final, norm_type='l2', degree_rise=2)
print("error l2:", error)
u_measured = Expression('sin(2*pi*(x[0]- cx*t))', degree=2, cx=Constant(0.1), t=0)
J = assemble((0.5 * (u_final - u_measured) ** 2) * dx)
dJdbeta0 = compute_gradient(J, [control])
h = Constant(0.1)
J_hat = ReducedFunctional(J, [control])
conv_rate = taylor_test(J_hat, beta_x, h) # check that J_hat is the derivative
print("convergence rate", conv_rate)
print("Gradient of dJdbeta0:", float(dJdbeta0[0]))
DJ[rk_ord].append(float(dJdbeta0[0]))
print(DJ)
def runge_kutta_convergence_oder():
rk_ord = 'RK2'
dg_ord = 2
m_sizes = [16, 32, 64, 128, 2 * 128, 4 * 128] #
errors = len(m_sizes) * [None]
for m, m_size in enumerate(m_sizes):
unit_mesh_, source_term, ic, beta, beta_x, gD = set_up_example(m_size)
u_final, t_final = forward(beta, gD, ic, source_term, unit_mesh_, h=1 / m_size, dg_order=dg_ord,
time_stepping=rk_ord)
errors[m] = errornorm(gD, u_final, norm_type='l2', degree_rise=2)
print(errors)
print(m_sizes)
convergence_rate = [errors[k - 1] / errors[k] for k in range(1, len(m_sizes))]
print(convergence_rate)
orders = [np.log(errors[k - 1] / errors[k]) / np.log((1 / m_sizes[k - 1]) / (1 / m_sizes[k])) for k in
range(1, len(m_sizes))]
print(orders)
if __name__ == '__main__':
runge_kutta_convergence_oder()
#dolfin_adjoint_taylor_test()
