Staggered approach to dolfin-adjoint problem

dolfin-adjoint
1

Hi all,
I’m trying to solve a problem that is linear in displacement, u and non linear in plastic strain with magnitude p and direction N. What I’d like to do is step up a form of the linear elasticity problem to begin with and solve for u at fixed p and N, then minimise a functional over p and N to enforce the plasticity rules. Since solving the variational elasticity problem is identical to minimising the elastic strain energy, I expect is that the result of this process should minimise a functional F = E_{strain} + \phi_{ajdoint}.

For now I’m just trying to get the first part to work, but something to do with the calculation of elastic strain is making things not work.

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

# Create mesh and define function space
length =1
height = 1
mesh = RectangleMesh(Point(0, 0), Point(length, height), 5, 5)

V_element = VectorElement('CG',mesh.ufl_cell(), 1)
T_element = TensorElement('DG',mesh.ufl_cell(), 0)
S_element = FiniteElement('DG',mesh.ufl_cell(), 0)                       
u_space = FunctionSpace(mesh, V_element)

u = Function( FunctionSpace(mesh, V_element))
p = Function( FunctionSpace(mesh, S_element))
N = Function( FunctionSpace(mesh, T_element))
u_ = TestFunction( FunctionSpace(mesh, V_element))

# Define elastic parameters
E = Constant(1.0e2)  # Young's modulus
nu = Constant(0.3)  # Poisson's ratio
Sy = Constant(0.3) # Yield stress
# Define the stress-strain relation
lambda_ = E * nu / ((1 + nu) * (1 - nu))
mu = E / (2 * (1 + nu))
I = Identity(2)

def eps(u): # elastic strain
    return sym(grad(u)) - N*p 
def sigma(u): #stress from elastic strain
    e = eps(u)
    return lambda_ * tr(e) * I + 2 * mu * e

# Define the variational problem
F = (inner(sigma(u), eps(u_)))* dx - inner(Constant((0,0)), u_) *dx

# Define Dirichlet boundary conditions
def left_boundary(x, on_boundary):
    return near(x[0], 0) and on_boundary
def right_boundary(x, on_boundary):
    return near(x[0], length) and on_boundary
def bottom_boundary(x, on_boundary):
    return near(x[1], 0) and on_boundary

bc_left = DirichletBC(u_space.sub(0), Constant(0.0), left_boundary)  # Fix left boundary
bc_bottom = DirichletBC(u_space.sub(1), Constant(0.0), bottom_boundary)  # Fix bottom boundary
bc_right = DirichletBC(u_space.sub(0), Constant(0.001), right_boundary)  # Fix right boundary
bcs = [bc_left, bc_right,bc_bottom]


# Solve variational problem
solve(F == 0, u, bcs)

Gives the following error:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "/root/Plasticity/SingleStep.py", line 67, in <module>
    solve(F == 0, u, bcs)
  File "/src/pyadjoint/fenics_adjoint/solving.py", line 51, in solve
    output = backend.solve(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 308, in _solve_varproblem
    problem = NonlinearVariationalProblem(eq.lhs, u, bcs, J,
  File "/usr/lib/python3/dist-packages/dolfin/fem/problem.py", line 157, in __init__
    F = Form(F, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
  File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 407, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 26, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 75, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).

I’m fairly sure it has something to do with the multiplication of N*p in the calculation of elastic strain, but I’m not sure why. All I’d like to do here is solve for u with fixed p and N, so that later on I can perform some other minimisation on p and N with a reduced functional. Is this possible or have I deeply misunderstood the dolfin-adjoint framework?

Kind regards,

Magnus

2

The problem here is that eps(u_) doesn’t have u_ in all arguments (i.e. the N*p term does not contain a test function. This does not make sense when you think about a variational form, as all terms is multiplied by a test function and then integrated over the integration domain.

3

Ah I see, I’m not sure how to remedy this though, as I can’t just add new test functions:

p_ = TestFunction( FunctionSpace(mesh, S_element))
N_ = TestFunction( FunctionSpace(mesh, T_element)) 
[...]
def eps(u,p, N):
    return sym(grad(u)) - N*p
[...]
F = (inner(sigma(u,p, N), eps(u_,p_, N_)))* dx - inner(Constant((0,0)), u_) *dx

As they’re not from the same space as u_:

Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_0
  v_0
  v_0

And I don’t think I can add an arbitrary u_ term to the definition of strain, as that would presumably invalidate the underlying equations of elasticity and mess with things later on. Do you have any suggestions @dokken ?

4

I don’t quite understand what you are trying to achieve here.
What is the strong formulation of the equation here? How does p and N look like in the PDE setting?

5

Sorry, I’m trying to use a variational approach to plasticity; the theory follows this thesis, and can be summed up by pages 53 and 148:
The state is defined by vector displacement u, plastic strain magnitude p and plastic strain direction N. (In future p will be \Delta p, tracking the plastic increment over multiple loading steps, but for now I’m just looking at a single step.) The whole method revolves around energy minimisation.

\epsilon =\frac{1}{2}( \nabla u + \nabla u^{T})
\epsilon_{elastic} = \epsilon - pN
\sigma = 2\mu_{lamé} \epsilon_{elastic} + \lambda_{lamé} tr(\epsilon_{elastic} )I

\psi_{elastic}(u,p,N) = \frac{1}{2} \sigma : \epsilon_{elastic}
So far so normal, but the way the yield criterion is enforced is via another potential term \phi(p) = \sigma_yp and some contraints on N, namely it must be traceless and N:N =\frac{2}{3}, which are enforced with lagrangian multipliers \lambda_1 and \lambda_2:
L = \lambda_1 tr(N) + \lambda_2 (N:N-2/3)
In the thesis this problem is said to be solved incrementally by minimising the functional J = (\psi_{elastic}(u,p,N) + \phi(p) +L) over [u,p,N,\lambda_1, \lambda_2] at each step (again, here I am just performing one step)
My approach so far has been to set up the elastic problem as the underlying PDE, effectively minimising \psi_{elastic} over u at a given p,N.
Then I intend to use \phi(p) +L as the functional to minimise over p,N,\lambda_1, \lambda_2 using dolfin-adjoint, with the bound p\ge 0. The result I hope to gain is a solution containing u, p, N that fall within the variables’s bounds and minimises J. I’ve been trying to get this to work in dolfinx but to no avail, which I think is to do with starting the problem on a boundary p=0, so the solvers couldn’t define a gradient and did not converge. I’ve moved to dolfin-adjoint as it looks like it handles bounded non-linear PDEs quite well.