# Rounding errors prevent progress

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 * 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},
# 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},
# method specific parameters:
"mem_lim": 100000,
"Hinit": "default",
})
file = File(
``````

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

I use version (2019.1.0).

Here is my boundaries,

``````from dolfin import *

def bcs_(mesh, V_vector, disp):

def boundary_bot(x, on_boundary):
tol = 1E-7
return on_boundary and near(x, 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.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", 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
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 * 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},
# method specific parameters:
"mem_lim": 100000,
"Hinit": "default",
})

sol = solver.solve()

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

opt_chi = optimization(chi)
``````

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

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