topopt_conical_main.py

import numpy as np
import nlopt
import resolution_functions
import optimization_functions
import matplotlib.pyplot as plt
import os
import meshio as mio
from config_topopt_conical_data import *

'''
© Simon Ans
File : topopt_conical_main.py

Optimize the permittivities on the mesh to maximize the reflection efficiencies
for different wavelengths on a given order.
'''

if __name__=="__main__":
    
    ## Prepare the computation ------------------------------------------------------------------------------------------------------------

    # Parallelization of the Finite Element Method
    os.system("export OPENBLAS_NUM_THREADS=%d"%num_threads_resolution)

    # Create the paramater file and the mesh
    resolution_functions.set_data()
    resolution_functions.generate_mesh()
    resolution_functions.make_computation_directories()
    
    ## Density field ----------------------------------------------------------------------------------------------------------------------

    msh_density_init = optimization_functions.init_density(res_folder_name + "template_density_table.pos", "rho_opt.pos", analytic_init_density)
    #msh_density_init = np.loadtxt("rho_opt.txt")

    if (Flag_test == 1):
        ## Test of the adjoint method -----------------------------------------------------------------------------------------------------
        if (Flag_design_metal == 1):
            optimization_functions.test_jacobian(msh_density_init, 7)
            #optimization_functions.test_jacobian_soloPoint(msh_density_init, 7, 40, 300)
        else :
            optimization_functions.test_jacobian(msh_density_init, 3)
            #optimization_functions.test_jacobian_soloPoint(msh_density_init, 3, 40, 350)
        resolution_functions.print_sum_energies()

    else :
        ## Optimization process -----------------------------------------------------------------------------------------------------------
        
        problem_size = np.size(msh_density_init)
        # Choose the MMA optimization (actually GCMMA from the 2002 paper of Svanberg)
        opt = nlopt.opt(nlopt.LD_MMA, problem_size)

        # Define the target
        opt.set_min_objective(optimization_functions.target)

        # Boundaries
        lb_array = low_bnd*np.ones(problem_size)
        ub_array = up_bnd *np.ones(problem_size)
        opt.set_lower_bounds(lb_array)
        opt.set_upper_bounds(ub_array)

        # Stopping criteria
        opt.set_ftol_abs(tol)    
        opt.set_maxeval(maxiter)

        if (Flag_binarize == 1):
            # Optimization with binarization
            for iter_bin in range(nb_binarizations+1):
                np.savetxt("iter_bin.txt", np.array([iter_bin]))
                print("######                Binarization iteration %d...            ######"%iter_bin)
                print("----------------------------------------------------------------------------------------------\n")
                rho_opt = opt.optimize(msh_density_init)
                msh_density_init = rho_opt.copy()
            optimization_functions.final_binarization(res_folder_name + "template_density_table.pos", "rho_opt.pos")
        else :
            np.savetxt("iter_bin.txt", np.array([0]))
            rho_opt = opt.optimize(msh_density_init)

        plot_msh_density(rho_opt, "rho_opt.pdf")
        plot_msh_density(np.loadtxt("jacobian.txt"), "jac_opt.pdf")
        
        # Print the optimization outcome
        sum_energies = resolution_functions.print_sum_energies()
        min_value = opt.last_optimum_value()
        successful_termination = opt.last_optimize_result()
        optimization_functions.feedback_result(min_value, successful_termination)