How to create a Functional graph in dolfin-adjoint?

dolfin-adjoint
1

In the example in the link below, I would like to make an evaluation, a graph, between the objective function value and the iteration number. I made a few attempts but without success. Can anybody help me?

How I transform, using the example of the link, the Functional J in a vector that has all the values of each iteration and can be used in a graph.

Here’s the link for the example:

http://www.dolfin-adjoint.org/en/latest/documentation/poisson-topology/poisson-topology.html

2

Dear,

I managed to generate a code that generates a curve between Functional J and the number of iterations. Now I would like some help to also insert the volume reduction in the graph. @dokken , could you help me?

I believe it could be the generation of a “vector” of the “volume_constraint”, but I couldn’t.

Follow the code below:

from fenics import *
from dolfin_adjoint import *
import ipopt
import matplotlib.pyplot as plt
import numpy as np

parameters["std_out_all_processes"] = False

Vol = Constant(0.4)  # volume bound on the control
p = Constant(5)  # power used in the solid isotropic material
eps = Constant(1.0e-3)  # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8)  # regularisation coefficient in functional

def k(a):
    """Solid isotropic material with penalisation (SIMP) conductivity
  rule, equation (11)."""
    return eps + (1 - eps) * a ** p

n = 250
mesh = UnitSquareMesh(n, n)
A = FunctionSpace(mesh, "CG", 1)  # function space for control
P = FunctionSpace(mesh, "CG", 1)  # function space for solution

class WestNorth(SubDomain):
    """The top and left boundary of the unitsquare, used to enforce the Dirichlet boundary condition."""

    def inside(self, x, on_boundary):
        return (x[0] == 0.0 or x[1] == 1.0) and on_boundary

bc = [DirichletBC(P, 0.0, WestNorth())]
f = interpolate(Constant(1.0e-2), P)  # the volume source term for the PDE

def forward(a):
    """Solve the forward problem for a given material distribution a(x)."""
    T = Function(P, name="Temperature")
    v = TestFunction(P)

    F = inner(grad(v), k(a) * grad(T)) * dx - f * v * dx
    solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
                                                              "maximum_iterations": 20}})

    return T

if __name__ == "__main__":
    a = interpolate(Vol, A)  # initial guess.
    T = forward(a)  # solve the forward problem once.
    
# Now we define the functional, compliance with a weak regularisation
# term on the gradient of the material

    J = assemble(f * T * dx + alpha * inner(grad(a), grad(a)) * dx)
    
    # Callback called at each iteration
    total_obj_list = []
    def eval_cb(j, a):
        total_obj_list.append(j)
    
    m = Control(a)
    Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)

    lb = 0.0
    ub = 1.0
    
    volume_constraint = UFLInequalityConstraint((Vol - a)*dx, m) # Volume Constraint
    problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=volume_constraint)
    parameters = {"acceptable_tol": 1.0e-3, "maximum_iterations": 100}
    solver = IPOPTSolver(problem, parameters=parameters)
    a_opt = solver.solve()
    
    print("The total objective throughout the optimization: ", total_obj_list)
    
    compliance = total_obj_list
    maxiter = len(compliance)
    plt.figure()
    plt.plot(np.arange(1, maxiter+1), compliance)
    plt.xlabel("Number of iterations")
    plt.ylabel("Compliance")
    plt.title("Convergence history", fontsize=16)
    plt.show();
    
    File("output/final_solution.pvd") << a_opt
    xdmf_filename = XDMFFile(MPI.comm_world, "output/final_solution.xdmf")
    xdmf_filename.write(a_opt)