Specific DOF of a Function as control variable

dolfin-adjoint
1

Hello,

Using dolfin_adjoint, is it possible to have, as control variables, some specific DOF of a function?
I know that we can have all of the DOF as control variables but what if I want to keep some of them fixed and only use the free ones as control variables?

Here is a minimal example showing what I would like to do:

from dolfin import *
from dolfin_adjoint import *                                                       
import numpy as np                                                                                                         

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

n = V.dim()
d = mesh.geometry().dim()                                 

# Recover the coordinates of the dof
dof_V_coordinates = V.tabulate_dof_coordinates()
dof_V_coordinates.resize((n, d))
dof_V_x = dof_V_coordinates[:, 0]                                                    
dof_V_y = dof_V_coordinates[:, 1]

# Some function as initialization                         
initial = Expression("sin(pi*x[0])*cos(pi*x[1])",degree=3)
u = interpolate(initial, V)

# find some indexes of the free and fixed dof
fixed_indexes = np.where(dof_V_x >= 0.5)[0]
free_indexes = np.where(dof_V_x < 0.5)[0]

fixed_value = 0.4
u.vector()[fixed_indexes] =  fixed_value


# Working example:
m = Control(u) #Control variables (all the dofs of the Function u)
J = assemble(inner(grad(u), grad(u)) * dx) # Some Functional 
Jhat = ReducedFunctional(J, m)
dJdm = Jhat.derivative()

# Non-working example:
m_1 = Control(u.vector()[free_indexes]) # Control variables (some dofs of the Function u)
Jhat_1 = ReducedFunctional(J, m_1)
dJdm_1 = Jhat_1.derivative()

I get: AttributeError: ‘numpy.ndarray’ object has no attribute ‘block_variable’ at the line where m_1 is declared. Does it mean that I cannot achieve this without using Custom functions? And if yes what are your suggestions?
To have it working, just comment the last three lines that are just there to better explain what I am trying to achieve.

Thank you for your time.

2

The simplest solution is to modify the gradient of the full problem, aka:

m = Control(u) #Control variables (all the dofs of the Function u)
J = assemble(inner(grad(u), grad(u)) * dx) # Some Functional
Jhat = ReducedFunctional(J, m)
dJdm = Jhat.derivative()
dJdm.vector()[fixed_indexes] = 0

However, this implies that you cannot use the built-in optimization methods.

To do this properly, you need to make an intermediate variable, which is the one that will be used in your expression. I have added a minimal working example of this strategy below:

from dolfin import *
from dolfin_adjoint import *
import numpy as np

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


# # Some function as initialization
initial = Expression("sin(pi*x[0])*cos(pi*x[1])",degree=3)
uh = interpolate(initial, V)
uh.rename("u","")
fixed_value = 0.4

class fixed_indices(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 0.5
mf = MeshFunction("size_t", mesh, 2, 0)
fixed_domain = fixed_indices()
fixed_domain.mark(mf, 1)
dx_c = Measure("dx", domain=mesh, subdomain_data=mf)
u, v = TrialFunction(V), TestFunction(V)
a = inner(u,v)*dx
l = inner(uh, v)*dx_c(0) + inner(Constant(fixed_value), v)*dx_c(1)

u = Function(V)
solve(a==l, u)

m = Control(uh) #Control variables
J = assemble(inner(grad(u), grad(u)) * dx) # Some Functional
Jhat = ReducedFunctional(J, m)
dJdm = Jhat.derivative()
print(dJdm.vector().get_local())

which yields

[ 0.2934843   0.60295222  0.07689859 -0.42947131 -0.1527415   0.90639199
  0.19458717  1.4761873  -0.52695545 -0.75811268  0.50790983 -0.03113176
 -0.99850619 -3.74725313 -0.39626439 -0.34956218  0.80677139 -0.0541857
 -0.38557608 -0.34352392  0.2649522   0.09599802 -1.75726921 -2.9235389
 -1.60314787  0.          0.          0.         -0.34876289 -0.13449154
  0.2601133   0.20486509 -0.92068632  0.33012212  0.0787908   0.31704804
 -1.19407037 -0.84793381  0.          0.          0.          0.
  0.          0.          0.         -0.63982149  0.56348174  0.37797151
 -0.02842044  1.38792627  0.9383983   0.19644893  0.          0.
  0.          0.          0.          0.          0.          0.
 -0.17723889  1.08410744  1.55911393  0.          0.          0.
  0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.           0.          0.]
1 Like
3

Thank you very much for your answer. I indeed wanted to use the built-in optimization methods.
I understand the method you proposed with the intermediate variable.
However, I was hoping for something less costly than solving a PDE each time I want to fix some DOFs. I plan on integrating it into some optimization based on a large number of degrees of freedom.

After quite some time, I think that I found a “more efficient” way to implement it with a custom function as introduced in dolfin_adjoint documentation.
Here is how I proceeded:

  • Custom function in a file/module named fix_dofs.py :
from fenics import *
from fenics_adjoint import *

def fix_dofs(func,fixed_indexes,fixed_values):             
    vec = func.vector()
    vec[fixed_indexes] = fixed_values
    return Function(func.function_space(), vec)
  • Overloading the function with the following module fix_dofs_overloaded.py:
from fenics import *
from fenics_adjoint import *
from pyadjoint import Block
from pyadjoint.overloaded_function import overload_function

from fix_dofs import fix_dofs


backend_fix_dofs = fix_dofs

class FixedDofsBlock(Block):
    def __init__(self, func, fixed_indexes, fixed_values, **kwargs):
        super(FixedDofsBlock, self).__init__() 
        self.kwargs = kwargs
        self.add_dependency(func)

        #variables that do not appear in the dependencies but are still used for the computation 
        self.fixed_indexes = fixed_indexes 
        self.fixed_values = fixed_values
        
    def __str__(self): 
        return 'FixedDofsBlock'

    def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None):
        adj_input = adj_inputs[0] 
        x = inputs[0].vector() 

        output = adj_input
        output[self.fixed_indexes] = 0.0 # the derivative of the fixed dofs is null
        return output
    
    def recompute_component(self, inputs, block_variable, idx, prepared):
        return backend_fix_dofs(inputs[0],self.fixed_indexes,self.fixed_values) 

fix_dofs = overload_function(fix_dofs, FixedDofsBlock)
  • The main script / MWE (Minimal Working Example):
from fenics import *
from fenics_adjoint import *
import numpy as np   
from numpy.random import rand

from fix_dofs_overloaded import fix_dofs

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'CG', 1)
uh = project(Expression("sin(pi*x[0])*cos(pi*x[1])",degree=3), V)

# Recover the coordinates of the dof
dof_V_coordinates = V.tabulate_dof_coordinates()
dof_V_x = dof_V_coordinates[:, 0]                                                    
dof_V_y = dof_V_coordinates[:, 1]
# find some indexes of the free and fixed dof
fixed_indexes = np.where(dof_V_x >= 0.5)[0]
free_indexes = np.where(dof_V_x < 0.5)[0]

fixed_value= 0.0
u = fix_dofs(uh,fixed_indexes,fixed_value)

J = assemble(inner(grad(u), grad(u)) * dx) # Some Functional 

m = Control(uh)
Jhat = ReducedFunctional(J, m)

h = Function(V)
h.vector()[:] = rand(h.vector().local_size())
taylor_test(Jhat, uh, h)