from fenics import * from fenics_adjoint import * from dolfin_adjoint import * import os import pandas as pd import numpy import time import matplotlib.pyplot as plt from shape_optimization import do_shape_opt import argparse def run(L_start,Lr,Ll,resolution,theta,material_dict,logdir,run_nr): C11,C12,C44 = material_dict['C11'],material_dict['C12'],material_dict['C44'] A = 2*C44/(C11-C12) pr = C12/C11 a1,a2,a3,a4 = material_dict['a1'],material_dict['a2'],material_dict['a3'],material_dict['a4'] rho0,rho1,rho2,rho3 = material_dict['rho0'],material_dict['rho1'],material_dict['rho2'],material_dict['rho3'] b_0,b_1,b_2,b_3 = material_dict['b_0'],material_dict['b_1'],material_dict['b_2'],material_dict['b_3'] delta = material_dict['delta'] start = time.time() set_working_tape(Tape()) E_end,at, ab, Delta, L_opt, chi_test, dpsi, u,verts, energy_dens,rotation = do_shape_opt(L_start,Lr,Ll,resolution,delta,theta,a1,a2,a3,a4,b_0,b_1,b_2,b_3,rho0,rho1,rho2,rho3) end = time.time() if not log_dir is None: datei = open(logdir+'/results.csv','a') datei.write('\n{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{}'.format(run_nr,Ll,L_start,Lr,E_end,at,ab,Delta,L_opt,end-start,theta,delta,resolution,verts,rho0,rho1,rho2,rho3,b_0,b_1,b_2,b_3,C11,C12,C44,A,pr)) datei.close() file = XDMFFile (log_dir+"/file_"+str(run_nr)+".xdmf") file.parameters["functions_share_mesh"] = True file.parameters ["rewrite_function_mesh"] = False energy_dens.rename("energy_dens","label") file.write(energy_dens,0) chi_test.rename("color","label") file.write(chi_test,0) dpsi.rename("psi","label") file.write(dpsi,0) u.rename("u","label") file.write(u,0) rotation.rename("rotation","label") file.write(rotation,0) file.close() if __name__ =="__main__": parser = argparse.ArgumentParser() parser.add_argument("--experiment", help="experiment to run, options are: single, materials, delta, anisotropy, poisson_ratio, theta. Defaults to single",default="single") parser.add_argument("--no_output", help="prevents output generation", action="store_true") args = parser.parse_args() experiment = str(args.experiment) no_output = args.no_output resolution = 150 theta = 0.25 CuAlNi = {"a1":11.41, "a2":-21.14, "a3":2.61, "a4":-10.94,"delta":0.0372,'rho0':0.5,'rho1':0.06,'rho2':0.2,'rho3':0.06,'b_0':1.,'b_1':45.9621,'b_2':0.192976,'b_3':13.1725,'C11':33.7,'C12':15.7,'C44':95} NiAl ={"a1":11.562724, "a2":-17.437087, "a3":10.062913, "a4":-9.375448,"delta":0.1,'rho0':0.5,'rho1':0.1,'rho2':0.2,'rho3':0.1,'b_0':1.,'b_1':45.1778,'b_2':3.95,'b_3':41.2111,'C11':115.5,'C12':45.5,'C44':110} Ybco = {"a1":7.19, "a2": -9.53, "a3":15.47, "a4":-1.88,"delta":0.01,'rho0':0.5,'rho1':0.1,'rho2':0.2,'rho3':0.1,'b_0':1.,'b_1':30.6768,'b_2':25.1363,'b_3':46.3434,'C11':209.,'C12':39.,'C44':100.} dicts = {'NiAl':NiAl,'CuAlNi':CuAlNi,'Ybco':Ybco} # create output directory if not no_output: if not os.path.isdir('output'):os.mkdir('output') log_dir = 'output/'+time.strftime("%b_%d_%H_%M_%S_", time.gmtime())+experiment if not os.path.isdir(log_dir):os.mkdir(log_dir) else: log_dir=None if not log_dir is None: datei = open(log_dir+'/results.csv','w') datei.write('{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{}'.format('run_nr','Ll','L_start','Lr','E','at','ab','Delta','L_opt','time','theta','delta','res','verts','rho0','rho1','rho2','rho3','b_0','b_1','b_2','b_3','C11','C12','C44','A','pr')) datei.close() if experiment == 'single': L_start = 6 Lr = 6 Ll = 2.5 theta = 0.25 material_dict = dicts['CuAlNi'] run(L_start,Lr,Ll,resolution,theta,material_dict,log_dir,0) elif experiment == 'materials': resolution = 250 theta = 0.25 Ll = 2.5 materials = ['NiAl','CuAlNi','Ybco'] L = 72.5 for run_nr, material_loc in enumerate(materials): material_dict = dicts[material_loc] L_start = 0.6*(1/material_dict['delta']) Lr = L-Ll-L_start run(L_start,Lr,Ll,resolution,theta,material_dict,log_dir,run_nr) elif experiment == 'delta': Ll = 2.5 theta = 0.25 material_dict = dicts['NiAl'] deltas = [0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.175,0.2] L = 22.5 for run_nr, delta_loc in enumerate(deltas): material_dict['delta'] = delta_loc L_start = 0.6*(1/delta_loc) Lr = L-Ll-L_start run(L_start,Lr,Ll,resolution,theta,material_dict,log_dir,run_nr) if not log_dir is None: fig = plt.plot(figsize=(10,10)) frame = pd.read_csv(log_dir+'/results.csv') frame.plot(x='delta', y='L_opt',logx=True,logy=True, marker='x') plt.savefig(log_dir+'/plot.pdf') plt.close() elif experiment == 'anisotropy': L_start = 6. Lr = 6 Ll = 2.5 theta = 0.25 resolution = 128 material_dict = dicts['NiAl'] C11 = material_dict['C11'] C12 = material_dict['C12'] C44 = material_dict['C44'] rho1=material_dict['rho1'] rho2=material_dict['rho2'] rho3=material_dict['rho3'] As = [0.7,0.8,0.9,1.,1.2,1.4,1.6,1.8,2.,2.25,2.5,2*C44/(C11-C12),4.,5.] C44s = [A*(C11-C12)*0.5 for A in As] for run_nr, C44_loc in enumerate(C44s): # solve equation (3.24) for different values of C44 b1,b2,b3 = numpy.dot(numpy.linalg.inv(numpy.identity(3)+numpy.array([[0,2*rho2,2*rho3],[rho1,rho2,2*rho3],[-rho1,-2*rho2,-3*rho3]])),numpy.array([C44_loc,(C11-C12)*0.5,C12])) material_dict['b_1'] = b1 material_dict['b_2'] = b2 material_dict['b_3'] = b3 material_dict['C44'] = C44_loc run(L_start,Lr,Ll,resolution,theta,material_dict,log_dir,run_nr) if not log_dir is None: fig = plt.plot(figsize=(10,10)) frame = pd.read_csv(log_dir+'/results.csv') frame.plot(x='A', y='L_opt', logx=True,logy=True, marker='x') plt.savefig(log_dir+'/plot.pdf') plt.close() elif experiment == 'poisson_ratio': L_start = 6 Lr = 6 Ll = 2.5 theta = 0.25 material_dict = dicts['NiAl'] C11 = material_dict['C11'] C12 = material_dict['C12'] C44 = material_dict['C44'] A = 2*C44/(C11-C12) rho1=material_dict['rho1'] rho2=material_dict['rho2'] rho3=material_dict['rho3'] prs = [0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,C12/C11,0.45,0.5] C12s = [pr*C11 for pr in prs] C44s = [A*(C11-C12_loc)*0.5 for C12_loc in C12s] for run_nr, C44_loc in enumerate(C44s): # solve equation (3.24) for different values of C44 b1,b2,b3 = numpy.dot(numpy.linalg.inv(numpy.identity(3)+numpy.array([[0,2*rho2,2*rho3],[rho1,rho2,2*rho3],[-rho1,-2*rho2,-3*rho3]])),numpy.array([C44_loc,(C11-C12s[run_nr])*0.5,C12s[run_nr]])) material_dict['b_1'] = b1 material_dict['b_2'] = b2 material_dict['b_3'] = b3 material_dict['C44'] = C44_loc material_dict['C12'] = C12s[run_nr] run(L_start,Lr,Ll,resolution,theta,material_dict,log_dir,run_nr) if not log_dir is None: fig = plt.plot(figsize=(10,10)) frame = pd.read_csv(log_dir+'/results.csv') frame.plot(x='pr', y='L_opt', marker='x') plt.savefig(log_dir+'/plot.pdf') plt.close() elif experiment == 'theta': L_start = 6 Lr = 6 Ll = 2.5 theta = 0.25 material_dict = dicts['NiAl'] thetas = [0.2,0.225,0.25,0.275,0.3,0.325,0.35,0.375,0.4,0.425,0.45,0.475,0.5] for run_nr, theta_loc in enumerate(thetas): run(L_start,Lr,Ll,resolution,theta_loc,material_dict,log_dir,run_nr) if not log_dir is None: fig = plt.plot(figsize=(10,10)) frame = pd.read_csv(log_dir+'/results.csv') frame.plot(x='theta', y='L_opt',logx=True,logy=True, marker='x') plt.savefig(log_dir+'/plot.pdf') plt.close() else: print ('experiment not found') sys.exit()