Optimizing Parameters in Initial Condition Expression

dolfin-adjoint
1

Hello all,

I’m trying to fit a simple reaction-diffusion model to data such that I identify a diffusion constant, degradation rate, and 3 parameters that define the initial conditions (concentration and position of a point source, modeled as a narrow gaussian).

When I minimize my error wrt my parameters, it seems to work for the diffusion constant and degradation term, but for the parameters that define my initial condition, I get:

WARNING:root:Adjoint value is None, is the functional independent of the control variable?

I assume this has to do with how I’m setting my initial conditions. I’ve also tried using “project” instead of “interpolate”. Can anyone give me advice?

Below is my full code, including the forward model and error estimation:

from fenics import *
from mshr import *
from fenics_adjoint import *
from collections import OrderedDict
import numpy as np
import matplotlib.pyplot as plt
from pyadjoint.overloaded_type import create_overloaded_object

#Define Geometry
R = 50.
domain = Circle(Point(0, 0), R)
mesh = generate_mesh(domain, 42)
mesh = create_overloaded_object(mesh)

T = 100.
# Define True Parameters:
D_A_true = Constant(float(np.random.random(1)))
g_true = Constant(.001*float(np.random.random(1)))

# BCs - PointSource + Neumann/NoFlux:
A_c0_true = Constant(float(np.random.randint(1,100)))
A_x_true = Constant(R)
A_y_true = Constant(R)
r = Constant(1.)

while np.sqrt((A_x_true(0))**2 + (A_y_true(0)**2))>R:
    A_x_true = Constant(R*(2*np.random.random()-1))
    A_y_true = Constant(R*(2*np.random.random()-1))
    
InputParams = [A_c0_true,A_x_true,A_y_true,D_A_true,g_true]


# Solve Forward Problem for True Solution
def forward(input_params,mesh_in):

    # Define Time:
    T = 100.
    num_steps = 10.
    dt = T/num_steps

    # Define Function Space
    V = FunctionSpace(mesh, "CG", 2)
    
    A_c0 = input_params[0]
    A_x = input_params[1]
    A_y = input_params[2]
    
    D_A = input_params[3]
    g = input_params[4]
    
    
    A = Function(V)
    ps = Expression('A_c0*exp(-(((x[0]-A_x)*(x[0]-A_x) + (x[1]-A_y)*(x[1]-A_y))/(r*r)))',A_c0=A_c0, A_x=A_x,A_y=A_y,r=r,degree=2)
    A_0 = interpolate(ps,V)
    assign(A,A_0)
    
    A_next = Function(V)
    v = TestFunction(V)

    F = ( inner((A_next - A)/dt,v) + D_A*inner(grad(A_next), grad(v)) + g*A_next*v)*dx

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

    As = {}
    
    t = 0.0
    while (t <= T):
        As[t] = A.copy()
        solve(F == 0, A_next,bc)
        A.assign(A_next)
        t += float(dt)
             
    return As



As_true = forward(InputParams,mesh)
#plot(As_true[T])
#plt.show()

# Define Test Parameters:
D_A = Constant(float(np.random.random(1)))
g = Constant(.001*float(np.random.random(1)))

# BCs - PointSource + Neumann/NoFlux:
A_c0 = Constant(float(np.random.randint(1,100)))
A_x = Constant(R)
A_y = Constant(R)
r = Constant(1.)

while np.sqrt((A_x(0))**2 + (A_y(0)**2))>R:
    A_x = Constant(R*(2*np.random.random()-1))
    A_y = Constant(R*(2*np.random.random()-1))

TestParams = [A_c0,A_x,A_y,D_A,g]

# Forward Function W/ Error Calculation
def forward_err(input_params,data,mesh_in):

# Define Time For Problem:
T = 100.
num_steps = 10.
dt = T/num_steps

V = FunctionSpace(mesh, "CG", 2)

# Initial Condition Parameters
A_c0 = input_params[0]
A_x = input_params[1]
A_y = input_params[2]

# System Parameters
D_A = input_params[3]
g = input_params[4]

# Setting Initial Conditions
A = Function(V)
ps = Expression('A_c0*exp(-(((x[0]-A_x)*(x[0]-A_x) + (x[1]-A_y)*(x[1]-A_y))/(r*r)))',A_c0=A_c0, A_x=A_x,A_y=A_y,r=r,degree=2)
A_0 = interpolate(ps,V)
assign(A,A_0)

A_next = Function(V)

v = TestFunction(V)

F = ( inner((A_next - A)/dt,v) + D_A*inner(grad(A_next), grad(v)) + g*A_next*v)*dx
bc = DirichletBC(V, 0.0, "on_boundary")

As = {}

t = 0.0
while (t <= T):
    As[t] = A.copy()
    solve(F == 0, A_next,bc)
    A.assign(A_next)
    t += float(dt)
         
j = assemble(((A - data))**2*dx)

    
return A, j

A_test,J = forward_err(TestParams,As_true[T],mesh)

m = [Control(c) for c in TestParams]
rf = ReducedFunctional(J, m)
# Minimize for each parameter
OptimizedParams = minimize(rf)
2

Please supply the complete code such that we can reproduce your error message.

3

Sorry! Just updated with the full code, including generation of the target pattern

4

First of all, overloading interpolation is not supported, and therefore you need to use project.
Second of all, if you use the Expression class, you need to supply the derivative manually.
See for instance this documented example.

However, a work-around is to use the SpatialCoordinate-function and projection.

    x = SpatialCoordinate(mesh)
    ps =A_c0*exp(-(((x[0]-A_x)*(x[0]-A_x) + (x[1]-A_y)*(x[1]-A_y))/(r*r)))
    A_0 = project(ps,V)
    assign(A,A_0)