I am currently trying to compute a consistent tangent for an elastic model using dolfin_adjoint and facing some issues that I cannot resolve.
As a background information: in continuum mechanics representative volume elements are frequently used to explore the microscopic behavior of a structure. By applying volume averages, the macroscopic behavior can be predicted. While going through the possibilities of accomplishing that task using FEniCS, I read about dolfin adjoint and its possibility to track the dependencies of a forward solution w.r.t specific parameters. I have only very little knowledge about adjoint equations but I thought: since one of the specific tasks can be the derivation of a macroscopic tangent - that means calculatiing the change of the averaged stress with the averaged deformation - dolfin_adjoint could be a useful tool. My research on how to use the tool made me think, that the approach presented in the code below should work:
# load required libraries
#########################
from dolfin import *
from dolfin_adjoint import *
from BoundaryConditions import PeriodicDomain, PointDomain
from PostProcessing import PostProcessingQuantity, PostProcessor, volumeAverage
import numpy as np
# Optimization options for the form compiler
############################################
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# relevant simulation parameters:
#################################
# domain information
domainBB=[[0, 0, 0],[1, 1, 0]] # domain bounding box
# incrementation and loading
eps11=Constant(0.5)
# time
nIncs=1
tMax=1 # maximimum simulation time
dt=1
# elastic parameters
E=2.1e11 # typical choice in our model
nu = 0.3
# Lame constants
lam = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))
# domain information (mesh, function spaces, ...)
#################################################
# import mesh
mesh = UnitSquareMesh(5, 5)
x=SpatialCoordinate(mesh)
dx=Measure("dx", domain = mesh, metadata={"quadrature_degree": 3})
# generate subdomains for periodicity constraints and fixed center point
periodicDomain=PeriodicDomain(mesh,domainBB,periodicAxes=[0,1],tolerance=1e-4) # periodicity constraints
# define function spaces
V = VectorElement("P", mesh.ufl_cell(), 2)
W = FunctionSpace(mesh,V,constrained_domain=periodicDomain)
# define relevant functions
u=Function(W)
w=TestFunction(W)
du=TrialFunction(W)
# define boundary conditions (fixed point)
BC=DirichletBC(W,Constant((0.,0.)),periodicDomain.fixedPoint,method="pointwise")
# required mathematical expressions
###################################
epsMacro = as_tensor([[eps11, 0.],[0., 0.]])
# Kinematics
############
dim = u.geometric_dimension()
I = Identity(dim)
epsilon = epsMacro + sym(nabla_grad(u))
# Kinetics
##########
sigma = 2*mu*epsilon + lam*tr(epsilon)*I # linear elastic model
# Define weak form
##################
WF = inner(sigma.T,grad(w))*dx
J = derivative(WF,u,du)
# Setup solver
##############
problem = NonlinearVariationalProblem(WF, u, BC, J=J,
form_compiler_parameters=ffc_options)
solver = NonlinearVariationalSolver(problem)
solverParams=solver.parameters
solverParams['newton_solver']['absolute_tolerance'] = 1e-2
solverParams['newton_solver']['relative_tolerance'] = 1E-6
solverParams['newton_solver']['maximum_iterations'] = 20
solverParams['newton_solver']['relaxation_parameter'] = 1.0
solverParams['newton_solver']['linear_solver'] = 'mumps'
# prepare post-processing
# ########################
# spaces
WVector = VectorFunctionSpace(mesh,"P",2)
WTensor = TensorFunctionSpace(mesh,"DG",1)
# get offset (coordinates of fixed point) for spatial coordinate to ensure "correct" boundary conditions
x0=as_vector(periodicDomain.fixedPoint.closestVertex.array()[0:dim])
# functions
postProcessor=PostProcessor()
postProcessor.add(dot(epsMacro,x-x0)+u, WVector, "Displacement")
postProcessor.add(epsilon, WTensor, "Strain")
postProcessor.add(sigma, WTensor, "Stress")
# Incremental solution
######################
t=0
nIncs=0
postProcessor.update(t)
solver.solve()
postProcessor.update(tMax)
# compute required volume averages
# epsHom=volumeAverage(epsilon,dx)
# sigHom=volumeAverage(sigma,dx)
#
# print(epsHom)
# print(sigHom)
m=Control(eps11)
V=assemble(1.*dx)
sig11=1/V*assemble(sigma[0,0]*dx)
C1111=compute_gradient(sig11, m)
This mimics a minimum working example of what I am trying to do. I know that the problem is linear: I just wanted to see if, when I plug in linear material behavior, I get out the same. This is also a reason for the application of a NonlinearVariationalSolver: my extension of this approach would require the framework to work with this kind of problem and solver.
However, whenever I try to execute this piece of code, I get the error:
File â/home/âŚpyadjoint/overloaded_type.pyâ, line 379, in _ad_annotate_block
block = self.block_class(*self._ad_args, **self._ad_kwargs)
TypeError: init() got an unexpected keyword argument âmethodâ
Does anybody know what happens here, or if I am completely wrong that what I am trying to do is manageable using this framework?
