UPDATED: Importing dolfin_adjoint causes change of behavior in forward solve

dolfin-adjoint
1

Hi all,
From my understanding, importing dolfin_adjoint should not modify the behavior of the forward model. However, upon importing dolfin_adjoint, ie:

from dolfin import *
from dolfin_adjoint import *

NonlinearVariationalSolver, which otherwise converges to the correct result suddenly fails. (The equations are from a nonlinear hyperelastic membrane problem, representing an air-filled membrane dam retaining water on one or both sides.) I have narrowed down the source of the problem to the following code (which updates a hydrostatic traction vector based on the current position of the membrane). This code is called inside a newton-raphson style load-increment loop, which repeatedly calls the NonlinearSolver.

...
# returns a Sum object describing the original position plus displacement vector (expression plus vector valued function)

position = membrane.get_position()   
pos = project(position, mem.V)
x_vals, _, z_vals = pos.split(True)
# x value corresponding to max z height
max_x = Constant(x_vals.vector()[np.argmax(z_vals.vector()[:])])
max_z = max(z_vals.vector()[:])
print(max_z)

At this point the value of max_z is no longer correct, which is the source of the solver failing (the traction becomes too large). Plotting position' andpos` confirms they are working correctly, therefore I am lead to believe dolfin_adjoint is somehow modifying the behavior of the above code snippet. Is this possible? Furthermore, the entire code and solver runs fine after importing dolfin_adjoint when this load is removed/replaced with a static (constant) traction.

The problem occurs with both 2019.1.0 and 2018.1.0 docker containers.
Any ideas are greatly appreciated.

I post the entire code for load below in case there is something obvious there that I missed:

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

class HydrostaticLoad(object):
    '''
    Hydrostatic (linearly varying) pressure field in the z direction.
    '''
    def __init__(self, membrane, rho=10, g=10, depth=.25, surround=False):

        self.membrane = membrane
        self.rho = Constant(rho)
        self.g = Constant(g)
        self.depth = Constant(depth)
        self.x_vec = Expression(('1', '0', '0'), degree=1)
        self.z_vec = Expression(('0', '0', '1'), degree=1)
        self.surround = surround

    def update(self):
        mem = self.membrane
        position = mem.get_position()

        x = dot(position, self.x_vec)
        z = dot(position, self.z_vec)

        pos = project(position, mem.V)
        x_vals, _, z_vals = pos.split(True)
        # x value corresponding to max z height
        max_x = Constant(x_vals.vector()[np.argmax(z_vals.vector()[:])])
        max_z = max(z_vals.vector()[:])
        print(np.argmax(z_vals.vector()[:]))
        print(max_z)  
        self.depth.assign(max_z)

        pressure = conditional(le(z,  self.depth), self.rho*self.g*(self.depth - z), 0.0)

        # return if "underwater" -- this is water on both sides
        if self.surround:
            pressure = project(pressure, mem.Vs)
            mem.p.assign(pressure)  # for visualization
            return pressure

        # Hack to get around situations where crown of membrane crosses over dry base
        pressure = conditional(le(x, max_x), pressure, 0)  # zero the pressure after the crown
        xi = SpatialCoordinate(mem.mesh)
        pressure = conditional(ge(xi[0], .9), 0, pressure)

        pressure = project(pressure, mem.Vs)
        return pressure
2

If you mean it changes the computations of the original code (and not that the replay is incorrect - which it likely will be because some of those operations are not currently annotated), then the only thing I see that dolfin-adjoint alters in the first code snippet posted is the pos.split(True) line.

The current overloaded split method assumes the deepcopy parameter is supplied through the keyword arguments. Thus, with dolfin-adjoint imported, the pos.split(True) will change to pos.split(). You could try pos.split(deepcopy=True) to check if this is the problem.

1 Like
3

Sebastian, thanks so much for your input. Indeed, passing in the deepcopy argument by keyword fixed the issue. Regarding the replay, could you please clarify which operations in the code above are not yet annotated? Thanks again!

4

So there are several functions that are not annotated here.
Function.split(deepcopy=True) is not annotated (I believe it isn’t very difficult to implement).
Constant(value) is not annotated (should be trivial to implement as it is already sort of implemented through the workaround c = Constant(0); c.assign(value)).
We do not support vector/numpy manipulations in the master branch. (I.e Function.vector() and Function.vector().get_local() / Function.vector()[:]). There is preliminary support for tracking to numpy arrays in the branch numpy_adjoint:
https://bitbucket.org/dolfin-adjoint/pyadjoint/branch/numpy_adjoint

Notably in this branch you should be able to get to work with slices of numpy arrays from Function vectors, but max and np.argmax are not implemented (yet). They might not be very difficult to implement - autograd could perhaps be used to implement the derivatives.

1 Like