Hello,
I’ve added my working code below, it’s a simple code but it’s not working correctly.
I want to do pde constraint optimization, I minimize the following functional, first solve my model for a prescribed displacement and then solve it for F_i, desired force, to minimize it. (I want my force-disp graph to pass through the desired force points after the optimization).
my functional:
I left the code running below, but the result never provides correctly. I also tried different solvers.
I would appreciate it if you could help me because I couldn’t get it to work somehow.
Thank you in advance!
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
import numpy as np
try:
import cyipopt
except ImportError:
print("""This example depends on IPOPT and pyipopt. \
When compiling IPOPT, make sure to link against HSL, as it \
is a necessity for practical problems.""")
raise
# Set log level
set_log_level(LogLevel.WARNING)
# 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}
l1 = 1e-1
l2 = 1e-9
normal = Constant((0, 1))
chi0 = 0.9
delta_t = 1e-1
E1, nu1 = 30.5, 0.4 # base material
E2, nu2 = 1e-10, 0.0 # void material
mu1, lmbda1 = Constant(E1 / (2 * (1 + nu1))), Constant(E1 * nu1 / ((1 + nu1) * (1 - 2 * nu1)))
mu2, lmbda2 = Constant(E2 / (2 * (1 + nu2))), Constant(E2 * nu2 / ((1 + nu2) * (1 - 2 * nu2)))
max_disp = [-0.1]
f_desired = [-0.11]
store_u = np.linspace(0, 1.0 * min(max_disp), 100)
displacement_by_step = store_u
eps = Constant(1e-6)
c1, c2 = 10, 1e-4
# Create mesh and define function space
n = 32
mesh = UnitSquareMesh(n, n)
W = FunctionSpace(mesh, "CG", 2)
V_vector = VectorFunctionSpace(mesh, "CG", 2)
W_tensor = TensorFunctionSpace(mesh, "CG", 2)
force_FS = VectorFunctionSpace(mesh, "Real", 0)
df = TestFunction(force_FS)
def bcs_(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
def H(a, l):
return 0.5 * (1 + a / (l + (a ** 2) ** 0.5 ))
# initial chi expression
chi0 = Constant(0.9)
chi = Function(W)
chi.interpolate(chi0)
file = File(
"./Results" + "/chi_initial.pvd");
file << chi
# Stored strain energy density (compressible neo-Hookean model)
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):
max_error = 0
count = 0
for i in max_disp:
u0 = Function(V_vector)
u = Function(V_vector)
step_size = np.linspace(0, i, 100, endpoint=True)
for j in step_size:
disp = Constant([0.0, j])
bcs, dss = bcs_(disp)
Ic, J_, v, F = kinematics(u)
# Stored strain energy density (compressible neo-Hookean model)
Psi = H((chi - 0.5), l2) * psi(mu1, Ic, J_, lmbda1) + (1 - H((chi - 0.5), l2)) * psi(mu2, Ic, J_, lmbda2)
# 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)
file2 = File("./Results" + "/stress-opt.pvd");
u.assign(u0)
solve(dPi == 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))
max_error += ((t - f_desired[count]) ** 2)
print("desired f", f_desired[count])
print("calculated f", t)
count += 1
J = c1 * max_error + assemble((((1/c2) * (chi * (1 - chi)) ** 2) + 0.5 * c2 * (inner(grad(chi), grad(chi)))) * dx)
print("max_error", max_error)
print("J", J)
controls = File("./Results" + "/control_iterations_guess.pvd")
chi_viz = Function(W, name="ControlVisualisation")
def eval_cb(j, chi):
chi_viz.assign(chi)
controls << chi_viz
control = Control(chi)
rf = ReducedFunctional(J, control, eval_cb_post=eval_cb)
problem = MinimizationProblem(rf) # , constraints=const
parameters_ipopt = {"acceptable_tol": 1.0e-10, "maximum_iterations": 100}
solver = IPOPTSolver(problem, parameters=parameters_ipopt) # , parameters=parameters_ipopt
f_opt = solver.solve()
file = File(
"./Results" + "/Xnew.pvd");
file << f_opt
plot(f_opt, title="f_opt")
plt.savefig('./Results' + '/Xnew' + '--' + str(i) + '.png')
print(f"Initial J: {rf(chi)}")
print(f"min opt {min(f_opt.vector().get_local())}, max opt {max(f_opt.vector().get_local())}")
print(f"Optimal J {rf(f_opt)}")
return f_opt
opt_chi = optimization(chi)
print("Optimization Done")
def forward(opt_chi, Traclist, u0):
u = Function(V_vector)
Ic, J_, v, F = kinematics(u)
# Stored strain energy density (compressible neo-Hookean model)
Psi = (H((opt_chi - 0.5), l2)) * psi(mu1, Ic, J_, lmbda1) + (1 - H((opt_chi - 0.5), l2)) * psi(mu2, Ic, J_, lmbda2)
# 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)
# Piola-Kirchhoff Stress
stress = diff(Psi, F)
# Solve variational problem
disp = Constant([0.0, 0.0])
bcs, dss = bcs_(disp)
recu = []
for d in displacement_by_step:
print(d)
disp.assign(Constant([0.0, d]))
u.assign(u0)
solve(dPi == 0, u, bcs)
u0.assign(u)
recu.append(u(0, 1))
file1 << u
stress_calculation = (project((stress), W_tensor))
traction = dot(stress, normal)
storedtraction = assemble((traction[1]) * dss(1))
Traclist.append(storedtraction)
return store_u, Traclist
Traclist = []
u0 = Function(V_vector)
store_u, Traclist = forward(H((opt_chi - 0.5), l2), Traclist, u0)
Traclist2 = []
u0 = Function(V_vector)
store_u2, Traclist2 = forward(H((chi - 0.5), l2), Traclist2, u0)
file = File(
"./Results" + "/chi_last.pvd");
file << project(H((opt_chi - 0.5), l2), W)
# Plot
K = -1
plt.figure()
displacement_by_step = [x * K for x in displacement_by_step]
Traclist = [y * K for y in Traclist]
plt.plot(displacement_by_step, Traclist, label="After- opt")
Traclist2 = [y * K for y in Traclist2]
plt.plot(displacement_by_step, Traclist2, label="Before- opt")
max_disp2 = [x * K for x in max_disp]
f_desired2 = [y * K for y in f_desired]
plt.plot(max_disp2, f_desired2, 'o', label="Desired Points")
plt.legend()
plt.xlabel("Displacement")
plt.ylabel("Force")
plt.title("Graph for the displacement" + str(max(max_disp)))
plt.savefig('./Results' + '/Graph.png')
I use a different chi expression from an image.