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


277
views
0
8 months ago by
danabl  
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)
    331         elif task_str.startswith(b'NEW_X'):
    332             # new iteration

/usr/lib/python2.7/dist-packages/scipy/optimize/lbfgsb.pyc in func_and_grad(x)
    277         def func_and_grad(x):
    278             f = fun(x, *args)
--> 279             g = jac(x, *args)
    280             return f, g
    281

/home/fenics/local/lib/python2.7/site-packages/dolfin_adjoint/optimization/optimization.pyc in <lambda>(m)
     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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »