Rounding errors prevent progress

dolfin-adjoint
1

Hi all,

I’m doing a pde constraint optimization. I have shared my minimal code below, I am getting a “Rounding errors prevent progress” error during optimization.

def optimization(chi, u0):
    for i in range(len(max_disp)):
        print("i:", i)
        disp = Constant([0.0, max_disp[i]])
        bcs, dss = bcs_(mesh, V_vector, disp)
        u = Function(V_vector)
        Ic, J_, v, F = kinematics(u)
        # Stored strain energy density (compressible neo-Hookean model)
        Psi = H(chi) * psi( Ic, J_, v, F)
        # Total potential energy
        Pi = Psi * dx
        # Compute first variation of Pi (directional derivative about u in the direction of v)
        dPi = derivative(Pi, u, v)
        step_size = np.linspace(0, max_disp[i], 100)
        solve(dPi == 0, u, bcs,
                  form_compiler_parameters=ffc_options)
        stress = diff(Psi, F)
        stress_calculation = (project(stress, W_tensor))
        t = dot(stress, normal)
        t = assemble(-t[1] * dss(1))
        error += (t - f_desired[i]) ** 2
    J =c1 * error + assemble( 0.5 * c2 * dot(grad(chi), grad(chi))) * dx)
    control = Control(chi)
    rf = ReducedFunctional(J, control)
    problem = MoolaOptimizationProblem(rf)
    f_moola = moola.DolfinPrimalVector(chi)
    moola.BFGS(problem, f_moola, options={'jtol': 0,
                                                   'gtol': 1e-9,
                                                   'Hinit': "default",
                                                   'maxiter': 200,
                                                   'mem_lim': 10})(rf, method="L-BFGS-B",
                        tol=1e-50, options={"jtol": 1e-50,
                                            "rjtol": 1e-50,
                                            "gtol": 1e-50,
                                            "rgtol": 1e-50,
                                            "maxiter": 200,
                                            "display": 2,
                                            "maxcorint": 1e+1000000,
                                            "line_search": "strong_wolfe",
                                            "line_search_options": {"ftol": 1e-50, "gtol": 1e-50,
                                                                    "xtol": 1e-1000,
                                                                    "start_stp": 1},
                                            "record": ("grad_norm", "objective"),
                                            # method specific parameters:
                                            "mem_lim": 100000,
                                            "Hinit": "default",
                                            })

    sol = solver.solve()

    f_opt = sol['control'].data
    return f_opt

Does using IPOPT at this point affect the result? Since there was a download problem, I changed my minimization as below, but this time I cannot get iterate enough.

f_opt = minimize(rf, method="L-BFGS-B",
                        tol=1e-50, options={"jtol": 1e-50,
                                            "rjtol": 1e-50,
                                            "gtol": 1e-50,
                                            "rgtol": 1e-50,
                                            "maxiter": 200,
                                            "display": 2,
                                            "maxcorint": 1e+1000000,
                                            "line_search": "strong_wolfe",
                                            "line_search_options": {"ftol": 1e-50, "gtol": 1e-50,
                                                                    "xtol": 1e-1000,
                                                                    "start_stp": 1},
                                            "record": ("grad_norm", "objective"),
                                            # method specific parameters:
                                            "mem_lim": 100000,
                                            "Hinit": "default",
                                            })
    file = File(
2

Please make sure that the code is executable by anyone copying it, by making sure that all imports and variables are defined.

Additionally, please supply what version of dolfin and dolfin-adjoint you are using.

3

I use version (2019.1.0).

Here is my boundaries,

from dolfin import *
from dolfin_adjoint import *

def bcs_(mesh, V_vector, disp):

    def boundary_bot(x, on_boundary):
        tol = 1E-7
        return on_boundary and near(x[1], 0.0, tol)
    
    
    bc_bot = DirichletBC(V_vector, [0.0, 0.0], boundary_bot)

    
    def boundary_top(x, on_boundary):
        tol = 1E-7
        return on_boundary and (near(x[1], 1.0, tol))
    
    bc_top = DirichletBC(V_vector, disp, boundary_top)
 
    bcs = [bc_bot, bc_top]
    boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundary_subdomains.set_all(0)
    AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)
    dss = ds(subdomain_data=boundary_subdomains)
    
    return bcs, dss

and here I shared the code again,

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

normal = Constant((0, 1))
Traclist = []

E1, nu1 = 26.0, 0.25  # base material
mu1, lmbda1 = Constant(E1 / (2 * (1 + nu1))), Constant(E1 * nu1 / ((1 + nu1) * (1 - 2 * nu1)))
max_disp = [-0.04, -0.09]
f_desired = [-0.04, -0.08]

n = 32
mesh = UnitSquareMesh(n, n, "crossed")
W0 = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 0)
V_vector = VectorFunctionSpace(mesh, "CG", 1)
W_tensor = TensorFunctionSpace(mesh, "CG", 1)
force_FS = VectorFunctionSpace(mesh, "Real", 0)
df = TestFunction(force_FS)
c1, c2 = 0.5, 1e-3
eps = 1e-6
chi = project(Expression("x[1]", name='Control', degree=3), W)
def H(a):
    return 0.5 * (1 + a / (1e-10 + (a **2 ) ** (1/2)))

def psi(mu, Ic, J_, lmbda):
    return (mu / 2) * (Ic - 3) - mu * ln(J_) + (lmbda / 2) * (ln(J_)) ** 2
 
 
def kinematics(u):
    v = TestFunction(V_vector)  # Test function
    d = u.geometric_dimension()
    I = Identity(d)  # Identity tensor
    F = variable(I + grad(u))  # Deformation gradient
    C = F.T * F  # Right Cauchy-Green tensor
    # Invariants of deformation tensors
    Ic = tr(C)
    J_ = det(F)
    return Ic, J_, v, F
def optimization(chi):
    for i in range(len(max_disp)):
        print("i:", i)
        disp = Constant([0.0, max_disp[i]])
        bcs, dss = bcs_(mesh, V_vector, disp)
        u = Function(V_vector)
        Ic, J_, v, F = kinematics(u)
        # Stored strain energy density (compressible neo-Hookean model)
        Psi =(eps - (1 - eps) * ((H((chi - 0.5), l2)) ** 3)) * psi(mu1, Ic, J_, lmbda1)
        # Total potential energy
        Pi = Psi * dx
        # Compute first variation of Pi (directional derivative about u in the direction of v)
        dPi = derivative(Pi, u, v)
        step_size = np.linspace(0, max_disp[i], 100)
        for j in step_size:
            disp = Constant([0.0, j])
            bcs, dss = bcs_(mesh, V_vector, disp)
            #udot = (u - u0) / delta_t
            #u.assign(u0)
            FF = dPi  # + nu * (inner(udot, v)) * dx
            solve(FF == 0, u, bcs,
                  form_compiler_parameters=ffc_options)
            # u0.assign(u)
            # Piola-Kirchhoff Stress
            stress = diff(Psi, F)
            stress_calculation = (project(stress, W_tensor))
            file2 << stress_calculation
            t = dot(stress, normal)
            t = assemble(-t[1] * dss(1))
        
        print("calculated force:", t)
        print("desired force:", f_desired[i])
        error += (t - f_desired[i]) ** 2
    J =c1 * error + assemble( 0.5 * c2 * dot(grad(chi), grad(chi))) * dx)
    control = Control(chi)
    rf = ReducedFunctional(J, control)
    problem = MoolaOptimizationProblem(rf)
    f_moola = moola.DolfinPrimalVector(chi)
    moola.BFGS(problem, f_moola, options={'jtol': 0,
                                                   'gtol': 1e-9,
                                                   'Hinit': "default",
                                                   'maxiter': 200,
                                                   'mem_lim': 10})(rf, method="L-BFGS-B",
                        tol=1e-50, options={"jtol": 1e-50,
                                            "rjtol": 1e-50,
                                            "gtol": 1e-50,
                                            "rgtol": 1e-50,
                                            "maxiter": 200,
                                            "display": 2,
                                            "maxcorint": 1e+1000000,
                                            "line_search": "strong_wolfe",
                                            "line_search_options": {"ftol": 1e-50, "gtol": 1e-50,
                                                                    "xtol": 1e-1000,
                                                                    "start_stp": 1},
                                            "record": ("grad_norm", "objective"),
                                            # method specific parameters:
                                            "mem_lim": 100000,
                                            "Hinit": "default",
                                            })

    sol = solver.solve()

    f_opt = sol['control'].data
    return f_opt


opt_chi = optimization(chi)
4

You will never any gradient satisfying this, as dolfin works with precision 1e-16.

6

Hi @dokken, do you have any suggestions to prevent this error? I am using default settings for jtol, rjtol, and gtol now.