General strategy for solving 1D energy minimization problem

Question

Let's say I have a non-linear energy functional E[u]. How do I solve for u that minimizes E over the whole mesh? For example

$$ E[u] = 2\int_0^W{E_d(z)u(z)^2dz} + \int_0^W{c(z)\left(f(z)-b(z)\int_0^z{u(v)^2dv}\right)^2dz} $$

This problem can be solved analytically resulting in a function that is constant below a certain $$ W_c < W $$ and zero above that.

Here is a (not so simple) example in dolfin

import dolfin
from dolfin import *
import ufl


# Optimization options for the form compiler
parameters['linear_algebra_backend'] = 'PETSc'
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_rule"] = 'auto'
ffc_options = {"optimize": True, \
           "eliminate_zeros": False, \
           "precompute_basis_const": True, \
           "precompute_ip_const": True
           }

class Primitive(ufl.Coefficient, dolfin.cpp.Expression):
    def __init__(self,f,mesh,element,constant= 0.0, **kwargs):
       dolfin.cpp.Expression.__init__(self)
       self._f = f
       self._mesh = mesh
       self._constant = constant
       subdomains = CellFunction('size_t',mesh)
       self._subdomains = subdomains
       dx = Measure('dx',domain=mesh)[subdomains]
       self._dx=dx
       self._F = Form(f*dx(1))

       ufl.Coefficient.__init__(self, element, count=self.id())

     def eval(self, values, x):

        inside_function = lambda y: y<=x[0]
        domain = AutoSubDomain(inside_function=inside_function)
        self._subdomains.set_all(0)
        domain.mark(self._subdomains,1)
        F = self._F
        values[0] = assemble(F) + self._constant

    def eval_cell(self, values, x, cell):
        self.eval(values,x)

    def ufl_evaluate(self,x, component, derivatives):
        values = np.zeros(1,dtype='d')
        self.eval(values,x)
        return values[0]


mesh = IntervalMesh(100,0,1e-6)
V = FunctionSpace(mesh, "CG", 1)

left = CompiledSubDomain("(std::abs(x[0]) < DOLFIN_EPS) && on_boundary")
right = CompiledSubDomain("(std::abs(x[0] - 1e-6) < DOLFIN_EPS) && on_boundary")

class LinearMismatchProfile(Expression):

    def __init__(self, G, W, W_b, f_0,**kwargs):
        self._W = W
        self._G = G
        self._W_b = W_b
        self._f_0 = f_0

    def eval(self, values, x):
         W = self._W
         W_b = self._W_b
         f_0 = self._f_0
         G = self._G

        if x[0] <= W-W_b:
            values[0] = G*x[0] + f_0
        else:
            values[0] = G*(W-W_b)

Ed = Constant(8.138931143006406e-09) #energy per unit length of dislocation
c = Constant(1e8) #Constant(92017082108.12033) #elastic constant for biaxial strain
b = Constant(2.161453136448972e-10) #part of Burger's vector responsible for relaxation (in-plane component)

W = 1e-6
W_b = 2.64e-7
G =  70000
f_0 = 0.0

f = LinearMismatchProfile(G=G,W=W,W_b=W_b,f_0=f_0,element=V.ufl_element())

u = Function(V) # dislocation density profile
v  = TestFunction(V) 
du = TrialFunction(V)

Pu = Primitive(f=u*u,mesh=mesh, constant=0.0, element=V.ufl_element())
epsilon = f - b*Pu #in-plane strain

E = 2*Ed*u*u*dx + c*(epsilon**2)*dx

F = derivative(E,u,v)
J = derivative(F,u,du) 


F = Form(F, form_compiler_parameters=ffc_options)
J = Form(J, form_compiler_parameters=ffc_options)

bcl = DirichletBC(V,sqrt(G/float(b)),left) 
bcr = DirichletBC(V,0.0,right)
bcs = []

# Set up the non-linear problem

class TersoffModelProblem(NonlinearProblem):
    def __init__(self, F, J, bcs=None, form_compiler_parameters=None, ident_zeros=False,
             exterior_facet_domains=None, interior_facet_domains=None):
        NonlinearProblem.__init__(self)
        self._F = F
        self._J = J
        self._bcs = bcs
        self._form_compiler_parameters = form_compiler_parameters


    def F(self, b, x):
        assemble(self._F, tensor=b,
             form_compiler_parameters=self._form_compiler_parameters,
             keep_diagonal=False)

        #Apply boundary conditions
        for bc in self._bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self._J, tensor=A,
             form_compiler_parameters=self._form_compiler_parameters,
             keep_diagonal=False)

        for bc in self._bcs:
            bc.apply(A)

        #if self.ident_zeros:
        #A.ident_zeros()

problem = TersoffModelProblem(F,J,bcs=bcs)


class TersoffModelSolver(NewtonSolver):

    def __init__(self, u, Pu, epsilon):
        NewtonSolver.__init__(self)
        self._epsilon = epsilon
        #self._epsilon_ufl = epsilon_ufl
        self._u = u
        self._Pu = Pu


#     def update_solution(self, x, dx):
#         NewtonSolver.update_solution(self,x,dx)
#         plot(self._u)
#         plot(self._Pu,mesh=mesh)

solver = TersoffModelSolver(u, Pu, epsilon)


solver_parameters={'nonlinear_solver': 'snes',
                      'newton_solver':
                      {"linear_solver": "lu",
                       "maximum_iterations": 100,
                       'relative_tolerance': 1e-30,
                       'absolute_tolerance': 1e-20,
                       'relaxation_parameter': 1.0,
                       },
                      'snes_solver':
                      {'sign': 'default',
                       'method': 'default',
                       'absolute_tolerance': 1e-12,
                       'solution_tolerance': 1e-12,
                       }
                     }

solver.parameters.update(solver_parameters['newton_solver'])

# Solve the problem
(iter, converged) = solver.solve(problem,u.vector())


# Plot and hold solution
plot(u, interactive = False)

This is similar to the strategy used for hyperelasticity demo. My attempt to use this strategy results in the trivial solution u=0. I'd like to know a general strategy to obtain the correct energy form.

Comments

You need to post a complete (but simple) example.

---

I just updated my question with a not so simple example because simple examples seem to converge. Any ideas?

Answer 1

For other readers, here's an example I wrote to show how to do global operators, as Chaffra is not the first to ask.

https://bitbucket.org/fenics-project/dolfin/issue/413/add-a-demo-of-global-operators-with-semi

Although in the code shown on this page there are some potential scaling problems, adding terms scaled by 1e8 and 1e-9 seems like asking for floating point precision trouble...

Comments

I tried to adpat your demo but it breaks when compiling the Jacobian. Please see below.

u = Function(V) # dislocation density profile
v  = TestFunction(V) 
du = TrialFunction(V)

F_form = u*dx(1) # TRYME
 # Pu = F(y, u) = \int^y{f(x)*dx(1)} where dx(1) depends on y and is defined by        mark_subdomain
Pu = FunctionalExpression(element=u.element())
Pu.init(F_form,update_functional)

# Compute partial derivative of Pu w.r.t. u.
# (NB! Note the difference between diff(u, u) == 1 and derivative(u, u, v) == v)
dFdu_form = diff(F_form, u)
dFdu = FunctionalExpression(element=u.element())
dFdu.init(dFdu_form, update_functional)

#Define in-plane strain coefficient
epsilon = f - b*Pu

#Define Energy density functional

 E = (2*Ed*u + c*epsilon**2)*dx

# Apply semi-automatic differentiation, by specifying the
# partial derivative of Pu w.r.t. u as another function dFdu:

#L1 = derivative(E, u, v, coefficient_derivatives={Pu: dFdu})


# Manually derived derivative of E
F = derivative(2*Ed*u*dx,u,v) + derivative(c*epsilon**2*dx, u, v, 
                                       coefficient_derivatives={Pu: dFdu})

J = derivative(F,u,du) 


F = Form(F, form_compiler_parameters=ffc_options)
J = Form(J, form_compiler_parameters=ffc_options)

Breaks here with

 /usr/local/lib/python2.7/dist-packages/ffc/analysis.pyc in _attach_integral_metadata(form_data, parameters)
259 
260     if all_equal(quad_schemes):
--> 261         scheme = quad_schemes[0]
262     else:
263         scheme = "canonical"

IndexError: list index out of range

---

You got a zero form after differentiation. The check was broken some time back, therefore the horrible error message. I can't spot the problem right now. Is it J or F that breaks?

---

It's J that breaks.

---

Pass the dict to all derivative calls. You may want to add other derivative relations to the dict as well.

---

OK that works!

---

Now I changed

F_form = u*u*dx(1)
F = derivative(2*Ed*u*u*dx,u,v) + derivative(c*epsilon**2*dx, u, v, 
                                   coefficient_derivatives={Pu: dFdu})

because the energy should always be positive, but NewtonSolver starts with zero residual so it gives me zero as the solution.

---

You say there is an analytical solution, try tracking down which functions and functionals compute the right values where you know what that is.

---

Yes I get a non-trivial solution when I use the analytical solution as my starting point. My problem now is that u must always be positive and so I am using the following hack to fullfil that constraints

class TersoffModelSolver(NewtonSolver):

    def __init__(self, u, Pu, epsilon, Pi):
        NewtonSolver.__init__(self)

        self._epsilon = epsilon
        self._u = u
        self._Pu = Pu
        self._Pi = Pi


    def update_solution(self, x, dx):

        i = (x-dx)<0.0
        dx[i] = 0.0

    NewtonSolver.update_solution(self,x,dx)

Although this gives me a solution close to my initial approximation it also lead to some oscillations on the solution probably because the Jacobian is not constrained. Any idea how to implement an inequality constraint in dolfin. I tried SNES but so far it gives me unphysical solutions when using sign='nonnegative'

---

I suggest creating a new question about inequality constraints, to attract people with such experience.