Defining Control variable with dolfin adjoint - no 'block_variable'


80
views
0
9 weeks ago by
Howdy,

Just trying to calculate sensitivity with dolfin adjoint:

solve(F == 0, uv, J=Jac)
u_r,u_i = split(uv)
obj_aug = ((u_r -cos_f)**2 + (u_i- sin_f)**2)*dx(2)
J =assemble(obj_aug)
dg = FunctionSpace(mesh,'DG',0)
w = Function(dg, name="Control")
gradient  = compute_gradient(J,Control(w))​

When I try to execute I get


self.block_variable = control.block_variable
AttributeError: 'Function' object has no attribute 'block_variable'
Abort trap: 6
​

Anyone know the quick fix to this? I'm not sure what this block variable is.
Community: FEniCS Project
Always include a minimum working example. If I can't run this, I can't help you. What's F?
written 9 weeks ago by Miguel  

Hi Miguel,

Thanks a lot for the comment. Cleaning up my code to post - I discovered that I had multiple files that imported dofin adjoint, but I imported * from my main file, which also had the import which apparently causes an issue. My entire code is here for folks that want an example for the future. Thank you Miguel!




from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import math
from random import randint
from numpy.random import rand
from petsc4py import PETSc
from mpl_toolkits.mplot3d import axes3d, Axes3D
from fenics_adjoint import *
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], -1.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], -1.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (-(5.0/8.0), (5.0/8.0))) and between(x[0], (-(5.0/8.0), (5.0/8.0))))
class Obstacletwo(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (-.5, .5)) and between(x[1], (.7, 1.0)))

class Obstaclethree(SubDomain):
    def inside(self, x, on_boundary):
        return (((between(x[0], (-1.0, -(5.0/8.0))) or between(x[0], ((5.0/8.0),1.0 )) ) and between(x[1], (-1.0, 1.0))) or (between(x[0], (-1.0, 1.0)) and (between(x[1], (-1.0, -(5.0/8.0)) or between(x[1], ((5.0/8.0), 1.0))))))

def helmholtz_state(w,x_size,y_size):
	left = Left()
	top = Top()
	right = Right()
	bottom = Bottom()
	obstacle = Obstacle()
	obstacletwo = Obstacletwo()
	obstaclethree = Obstaclethree()

	# Create mesh and define function space
	p1 = Point(-1.0,-1.0)
	p2 = Point(1.0,1.0)
	mesh = RectangleMesh(p1,p2,x_size,y_size)
	element = VectorElement('P', triangle, 1, dim=5)
	V = FunctionSpace(mesh, element)
	W= FunctionSpace(mesh, "CG", 1)
	K = FunctionSpace(mesh,"DG",0)

	# Initialize mesh function for interior domains
	domains = MeshFunction("size_t", mesh,mesh.topology().dim())
	domains.set_all(0)
	obstacle.mark(domains, 1)
	obstacletwo.mark(domains, 2)
	obstaclethree.mark(domains, 3)



	# Initialize mesh function for boundary domains
	boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
	boundaries.set_all(0)
	left.mark(boundaries, 1)
	top.mark(boundaries, 2)
	right.mark(boundaries, 3)
	bottom.mark(boundaries, 4)

	dx = Measure("dx")(subdomain_data=domains)
	ds = Measure("ds")(subdomain_data=boundaries)



	# Define variational problem
	
	cos_f = Expression("cos(6*pi*x[1])",degree=2)
	sin_f = Expression("sin(6*pi*x[1])",degree=2)


	element1 = VectorElement('CG', mesh.ufl_cell(),2)
	V1 = FunctionSpace(mesh, element1)
	dg = FunctionSpace(mesh,'DG',0)
	wp = FunctionSpace(mesh,'CG',1)
	W = FunctionSpace(mesh,'CG',1)

	w= Function(dg)

	
	
	
	
	
	uo_r = Expression("36*pi*pi*cos(6*pi*x[1])",degree=2)
	uo_i = Expression("36*pi*pi*sin(6*pi*x[1])",degree=2)
	ko= Constant(6*pi)
	kosquared= Constant(36*pi*pi)
	q= Constant(.75)
	uv=Function(V1)
	u_r,u_i = split(uv)
	v_r,v_i = TestFunction(V1)
	F = inner(grad(u_r),grad(v_r))*dx + inner(grad(u_i),grad(v_i))*dx + ko*u_i*v_r*ds(1) -ko*u_r*v_i*ds(1)+ ko*u_i*v_r*ds(2) -ko*u_r*v_i*ds(2)+ ko*u_i*v_r*ds(3) -ko*u_r*v_i*ds(3)+ ko*u_i*v_r*ds(4) -ko*u_r*v_i*ds(4)  - kosquared*u_r*v_r*dx -kosquared*q*w*u_r*v_r*dx(1) -kosquared*u_i*v_i*dx -kosquared*q*w*u_i*v_i*dx(1)  - uo_r*q*w*v_r*dx(1) - uo_i*q*w*v_i*dx(1)  

	Jac = derivative(F, uv, TrialFunction(V1))
	
	
	
	# # break
	solve(F == 0, uv, J=Jac)
	u_r,u_i = split(uv)
	obj_aug = ((u_r -cos_f)**2 + (u_i- sin_f)**2)*dx(2)
	J =assemble(obj_aug)
	print(type(w))
	gradient  = compute_gradient(J,Control(w))
	
	# dolfin.plot(gradient)
	# plt.show()




for i in range(5,6):
	helmholtz_state(None,2**i,2**i)
written 9 weeks ago by Ryan Vogt  

1 Answer


0
9 weeks ago by
Summary: I'm a moron that doesn't know how to import correctly. If someone wants an example here it is.

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import math
from random import randint
from numpy.random import rand
from petsc4py import PETSc
from mpl_toolkits.mplot3d import axes3d, Axes3D
from fenics_adjoint import *
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], -1.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], -1.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (-(5.0/8.0), (5.0/8.0))) and between(x[0], (-(5.0/8.0), (5.0/8.0))))
class Obstacletwo(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (-.5, .5)) and between(x[1], (.7, 1.0)))

class Obstaclethree(SubDomain):
    def inside(self, x, on_boundary):
        return (((between(x[0], (-1.0, -(5.0/8.0))) or between(x[0], ((5.0/8.0),1.0 )) ) and between(x[1], (-1.0, 1.0))) or (between(x[0], (-1.0, 1.0)) and (between(x[1], (-1.0, -(5.0/8.0)) or between(x[1], ((5.0/8.0), 1.0))))))

def helmholtz_state(w,x_size,y_size):
	left = Left()
	top = Top()
	right = Right()
	bottom = Bottom()
	obstacle = Obstacle()
	obstacletwo = Obstacletwo()
	obstaclethree = Obstaclethree()

	# Create mesh and define function space
	p1 = Point(-1.0,-1.0)
	p2 = Point(1.0,1.0)
	mesh = RectangleMesh(p1,p2,x_size,y_size)
	element = VectorElement('P', triangle, 1, dim=5)
	V = FunctionSpace(mesh, element)
	W= FunctionSpace(mesh, "CG", 1)
	K = FunctionSpace(mesh,"DG",0)

	# Initialize mesh function for interior domains
	domains = MeshFunction("size_t", mesh,mesh.topology().dim())
	domains.set_all(0)
	obstacle.mark(domains, 1)
	obstacletwo.mark(domains, 2)
	obstaclethree.mark(domains, 3)



	# Initialize mesh function for boundary domains
	boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
	boundaries.set_all(0)
	left.mark(boundaries, 1)
	top.mark(boundaries, 2)
	right.mark(boundaries, 3)
	bottom.mark(boundaries, 4)

	dx = Measure("dx")(subdomain_data=domains)
	ds = Measure("ds")(subdomain_data=boundaries)



	# Define variational problem
	
	cos_f = Expression("cos(6*pi*x[1])",degree=2)
	sin_f = Expression("sin(6*pi*x[1])",degree=2)


	element1 = VectorElement('CG', mesh.ufl_cell(),2)
	V1 = FunctionSpace(mesh, element1)
	dg = FunctionSpace(mesh,'DG',0)
	wp = FunctionSpace(mesh,'CG',1)
	W = FunctionSpace(mesh,'CG',1)

	w= Function(dg)

	
	
	
	
	
	uo_r = Expression("36*pi*pi*cos(6*pi*x[1])",degree=2)
	uo_i = Expression("36*pi*pi*sin(6*pi*x[1])",degree=2)
	ko= Constant(6*pi)
	kosquared= Constant(36*pi*pi)
	q= Constant(.75)
	uv=Function(V1)
	u_r,u_i = split(uv)
	v_r,v_i = TestFunction(V1)
	F = inner(grad(u_r),grad(v_r))*dx + inner(grad(u_i),grad(v_i))*dx + ko*u_i*v_r*ds(1) -ko*u_r*v_i*ds(1)+ ko*u_i*v_r*ds(2) -ko*u_r*v_i*ds(2)+ ko*u_i*v_r*ds(3) -ko*u_r*v_i*ds(3)+ ko*u_i*v_r*ds(4) -ko*u_r*v_i*ds(4)  - kosquared*u_r*v_r*dx -kosquared*q*w*u_r*v_r*dx(1) -kosquared*u_i*v_i*dx -kosquared*q*w*u_i*v_i*dx(1)  - uo_r*q*w*v_r*dx(1) - uo_i*q*w*v_i*dx(1)  

	Jac = derivative(F, uv, TrialFunction(V1))
	
	
	
	# # break
	solve(F == 0, uv, J=Jac)
	u_r,u_i = split(uv)
	obj_aug = ((u_r -cos_f)**2 + (u_i- sin_f)**2)*dx(2)
	J =assemble(obj_aug)
	print(type(w))
	gradient  = compute_gradient(J,Control(w))
	
	# dolfin.plot(gradient)
	# plt.show()




for i in range(5,6):
	helmholtz_state(None,2**i,2**i)​

Please login to add an answer/comment or follow this question.

Similar posts:
Search »