Derivative of Augmented Lagrangian term

General
1

Hi, friends.

I’m struggling to implement the derivative of a term in my Augmented Lagrangian formulation.

This is the term representing the summation of the constraints and the sensitivities:

derivative
derivative929×239 14.1 KB

I’m verifying the sensitivities using the finite difference method, but the results are incorrect. I can’t figure out what I’m doing wrong. I’m using the following code to perform both calculations:

    def constraint(self):
        """
        Stress constraint.
        """

        # Short notation
        d_phi = self.par.delta_phi_test
        dx = self.par.dx
        H = self.par.H
        lb = self.par.lambda_2
        p = self.par.p
        r = self.par.r
        s_vm = self.par.von_mises
        s_y = self.par.sigma_y

        # Define condition
        G_bar = lb/r + H**p * s_vm / s_y - 1
        condition = ufl.conditional(ufl.gt(G_bar, 0), G_bar, 0)

        # Calculate the constraint integral at each node
        constraint_integral = assemble(r/2 * condition**2 * d_phi * dx)

        # Shape functions
        shape = assemble(d_phi * dx)

        # Calculate the constraint in each node
        constraint_nodal = constraint_integral.get_local() / shape.get_local()

        return constraint_nodal.sum()

   def Dconstraint_dphi(self):
        """
        Derivative of the stress constraint
        """

        # Create the indices
        i, j, _ = ufl.indices(3)

        # Short notation
        d_phi = self.par.delta_phi_test
        DH = self.par.DH
        dx = self.par.dx
        H = self.par.H
        lb = self.par.lambda_2
        p = self.par.p
        P = self.par.P
        r = self.par.r
        s = self.par.sigma
        s_vm = self.par.von_mises
        s_y = self.par.sigma_y

        # Auxiliary matrix V0
        V = as_tensor([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, 3]])

        # Define the condition
        G_bar = lb/r + H**p * s_vm / s_y - 1
        condition = ufl.conditional(ufl.gt(G_bar, 0), G_bar, 0)

        # Summation
        sum_ = (r * condition * 1./s_y * 1./2. * (s[i] * V[i, j] * s[j])**(
            -1./2.) * (V[i, j] * s[j] + V[j, i] * s[j]) * p * H**(p-1.) * s[i]
            * DH)
        sum_per_node = r * assemble(sum_ * dx).get_local(
        ) / assemble(d_phi * dx).get_local()
        self.par.sum_per_node.vector()[:] = sum_per_node[:]
        self.par.sum_per_node.vector().apply('insert')

        # Calculate the derivative of G in relation to phi
        Dconstraint = assemble(self.par.sum_per_node * d_phi * dx)

        # Calculate the derivative of G in relation to phi
        return Dconstraint.get_local()

I feel confident that the math in the paper is correct, but I don’t know how to properly implement it. My code for verifying sensitivities is working correctly, as the other sensitivities match.

I would really appreciate any kind of help.

2

As the code is very far away from reproducible, it is very hard to give you any guidance (as this is cleary just one function taken out of a bigger class).

This means that there are tons of undefined variables, that one can only guess what is.

Secondly, without adding the reference to the paper that you are taking the mathematics from, it is quite hard to get the context of all the terms.

3

Thanks for the answer, Dokken.

Actually, I calculated it by hand. Despite it being a well-known approach, I’m not basing it on any paper, so some equations can’t be found in any publication (the smoothed Heaviside, for example, is a sigmoid that I formulated).

This is the first time I’m asking for help here, and I’m really unsure of the best approach to take.

How can I prepare my question properly? I really appreciate any kind of help.

4

One would preferably need a minimal reproducible example.
See the guidelines: