### How to use dolfin-adjoint for optimization on multiple subdomains?

277
views
0
8 months ago by
Hi

I would like to use dolfin-adjoint for data assimilation problems in which model parameters are expected to assume different constant values on different subdomains. As an example, assume a diffusion problem in which the domain is divided into two subdomains O_1, O_2 with diffusion parameters D_1 and D_2 respectively.

I managed to set up forward models with multiple subdomains, as well as (multi-parametric) dolphin-adjoint optimization problems on a single domain. However, I am struggling to get multi-domain optimization working.

The following MWE specifies a 2D diffusion problem with 2 different diffusion coefficients, one on each half of the domain. Model parameters are assigned to the different subdomains by means of an expression class.
from dolfin import *
from dolfin_adjoint import *

#----- EXPRESSION CLASS FOR PARAMETERS ON SUBDOMAINS
class DiscontinuousScalar(Expression):
def __init__(self, cell_function, scalars,**kwargs):
self.cell_function = cell_function
self.coeffs        = scalars

def eval_cell(self, values, x, cell):
subdomain_id = self.cell_function[cell.index]
local_coeff  = self.coeffs[subdomain_id]
local_coeff.eval_cell(values, x, cell)

#----- MESH
nx = ny = 20
mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)

#----- MESH SUBDOMAINS
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0]-x[1] <= 0.0 + FENICS_EPS

class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0]-x[1] >= 0.0 - FENICS_EPS

materials = CellFunction('size_t', mesh)
materials.set_all(0)
subdomains = [Omega_0(), Omega_1()]
for m, subdomain in enumerate(subdomains,0):
subdomain.mark(materials, m)

#----- PROBLEM DEFINITION
V = FunctionSpace(mesh, "P", 1)
u_0_expression = Expression('exp(-a*pow(x[0], 2) - a*pow(x[1], 2))', degree=2, a=2)
u_0 =  interpolate(u_0_expression, V)
bc = DirichletBC(V, 0.0, "on_boundary")
time_step = 0.01
n_steps   = 100

def solve_diffusion(params):
D1 = params[0]
D2 = params[1]
D = DiscontinuousScalar(materials, [D1, D2], degree=1)
dt  = Constant(time_step)
u_n = u_0.copy(deepcopy=True)
u   = TrialFunction(V)
v   = TestFunction(V)
F = u*v*dx + dt*D*dot(grad(u), grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)
u = Function(V)
t = 0
for n in range(n_steps):
t += float(dt)
solve(a == L, u, bc)
u_n.assign(u)
return u

#----- SOLUTION FORWARD PROBLEM
D1 = Constant(1.2)
D2 = Constant(0.5)
u_target = solve_diffusion([D1,D2])

#----- OPTIMISATION
D1 = Constant(0.0001)
D2 = Constant(0.1)
u = solve_diffusion([D1,D2])
J = Functional(inner(u-u_target, u-u_target)*dx*dt[FINISH_TIME])
reduced_functional = ReducedFunctional(J, [Control(D1),Control(D2)])
m_opt = minimize(reduced_functional, method = "L-BFGS-B",
tol=2e-08, options = {"disp": True})
u_opt = solve_diffusion(m_opt)​

The forward model runs fine; and also parameter optimisation works on a simplified model (single domain). In this MWE, however, optimization throws the following error (see full error log at bottom of post):
TypeError: 'NoneType' object is not iterable

I assume this is somehow related to the way how parameters D_1 and D_2 are assigned to the different subdomains...?

What is the best approach for
1. specifying model parameters on subdomains (ideally read from xml mesh files), in a way that is
2. compatible with dolfin-adjoint optimization?

Thanks for you help!
best,
Daniel

Full error log:
Using dolfin-version 2017.1.0 from the latest dolfin-adjoint docker image (quay.io/dolfinadjoint/dolfin-adjoint)
Machine precision = 2.220D-16
N =            2     M =           10
This problem is unconstrained.

At X0         0 variables are exactly at the bounds
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/home/fenics/fenics-dev/diffusion_multipleDomain_optimisation_mwe.py in <module>()
77 reduced_functional = ReducedFunctional(J, [Control(D1),Control(D2)])
78 m_opt = minimize(reduced_functional, method = "L-BFGS-B",
---> 79                  tol=2e-08, options = {"disp": True})
80 u_opt = solve_diffusion(m_opt)
81

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/misc.pyc in func_wrapper(*args, **kwargs)
31    def func_wrapper(*args, **kwargs):
32        with annotations(False):
---> 33            res = func(*args, **kwargs)
34        return res
35

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/optimization/optimization.pyc in minimize(rf, method, scale, **kwargs)
221         kwargs["method"] = method
222
--> 223     opt = algorithm(rf_np, **kwargs)
224
225     if len(opt) == 1:

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/optimization/optimization.pyc in minimize_scipy_generic(rf_np, method, bounds, **kwargs)
126         res = scipy_minimize(J, m_global, method=method, bounds=bounds, **kwargs)
127     else:
--> 128         res = scipy_minimize(J, m_global, method=method, **kwargs)
129
130     rf_np.set_controls(np.array(res["x"]))

/usr/lib/python2.7/dist-packages/scipy/optimize/_minimize.pyc in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
445     elif meth == 'l-bfgs-b':
446         return _minimize_lbfgsb(fun, x0, args, jac, bounds,
--> 447                                 callback=callback, **options)
448     elif meth == 'tnc':
449         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,

/usr/lib/python2.7/dist-packages/scipy/optimize/lbfgsb.pyc in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)
328                 # minimization routine wants f and g at the current x
329                 # Overwrite f and g:
--> 330                 f, g = func_and_grad(x)
332             # new iteration

278             f = fun(x, *args)
--> 279             g = jac(x, *args)
280             return f, g
281

60     m_global = rf_np.obj_to_array(m)
61     J = rf_np.__call__
---> 62     dJ = lambda m: rf_np.derivative(m, forget=forget, project=project)
63     H = rf_np.hessian
64

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/reduced_functional_numpy.pyc in derivative(self, m_array, forget, project)
74             self(m_array)
75
---> 76         dJdm = self.__base_derivative__(forget=forget, project=project)
77
78         if project:

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/reduced_functional.pyc in derivative(self, forget, project)
260
261         # Apply the scaling factor
--> 262         scaled_dfunc_value = [utils.scale(df, self.scale) for df in list(dfunc_value)]
263
264         # Call callback

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/utils.pyc in scale(obj, factor)
34         except TypeError:
35             # Lists of dolfin objects?
---> 36             scaled_obj = [scale(o, factor) for o in obj]
37
38     return scaled_obj

TypeError: 'NoneType' object is not iterable​

Community: FEniCS Project