# Computation of conistent tangent using dolfin_adjoint

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 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

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)

# 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)

# Kinetics
##########
sigma = 2*mu*epsilon + lam*tr(epsilon)*I                                            # linear elastic model

# Define weak form
##################
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()

# 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)
``````

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:

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?

1 Like

Just as additional information: I am using FEniCS 2019.1.0 that I installed using anaconda and installed dolfin_adjoint 2019.1.0 via PiP.

Hi Flifff, I am currently facing the same problem and wanted to know if you by change have found a solution to this issue?

As I use DG to solve my PDE I must use âgeometricâ boundary conditions

`Preformatted textbcp_outflow = DirichletBC(V_dg, u_out, outflow, method='geometric')`

I get the same error if I run the code using

``````from fenics import *
``````

Does anyone know how the method `DirichletBC` is overloaded in Dolfin-Adjoint?

Hi RebeccaE,

unfortunately I gave up on this topic, for now. I did not manage to make it work and have continued my simulations without using dolfin_adjoint. Sorry, I hope youâll find a solution.

1 Like

Itâs typically iladvised to use strong imposition of boundary conditions with DG FE formulations. Consider weakly imposing those Dirichlet BCs. This will also remove any requirement for you to alter the underlying matrix with a `DirichletBC` object and the complications you see with `dolfin_adjoint`. See dolfin_dg for example if you have a complicated FE form.

2 Likes

The DirichletBC method error should not occur on the dolfin-adjoint master branch, which you can find here: https://github.com/dolfin-adjoint/pyadjoint

You can install it directly by using

``````pip install git+https://github.com/dolfin-adjoint/pyadjoint.git
``````
3 Likes

Thank you all for your help!

I will test the two suggestions and write what worked

Using the version of the master solved the error for me. Thank you for the good tip!

1 Like

Hi nate,

thank you for the link to the dolfin dg demos

Iâm not quite sure what the main problem is with strongly imposed boundary conditions for DG FE. Could you shortly explain?

How one would impose the DirecletBC weakly e.g. for the demo example demo_dg-advection-diffusion? (or any other simple demo example)