From 7bdf04b9e559f6f0cd6ea960fd715e714fa21285 Mon Sep 17 00:00:00 2001 From: Guillaume Demesy <guillaume.demesy@fresnel.fr> Date: Tue, 14 Feb 2023 21:05:44 +0100 Subject: [PATCH] new loop over theta --- DiffractionGratings/grating3D.geo | 5 ++ DiffractionGratings/grating3D_data_pillar.geo | 10 +-- .../grating3D_parallel_Mmatrix.sh | 65 ++++++++------- .../grating3D_postplot_Mmatrix.py | 79 ++++++++++--------- 4 files changed, 90 insertions(+), 69 deletions(-) diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo index ea1907e..8241901 100644 --- a/DiffractionGratings/grating3D.geo +++ b/DiffractionGratings/grating3D.geo @@ -213,6 +213,11 @@ Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz}; EndIf + If (tag_geom==8) + Cylinder(9) = {0,0,hh_L_3 , 0,0,rz , rx}; + Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, 1 } } { Volume{9}; } + EndIf + Coherence; Call SetPBCs; diff --git a/DiffractionGratings/grating3D_data_pillar.geo b/DiffractionGratings/grating3D_data_pillar.geo index 71b45aa..8871b78 100644 --- a/DiffractionGratings/grating3D_data_pillar.geo +++ b/DiffractionGratings/grating3D_data_pillar.geo @@ -4,7 +4,7 @@ // solve wall time pastix : 966.046s - R00=0.2432916553175064 // beware : very near a Rayleigh anomaly -nm = 1000; +nm = 1; pp1 = "1Incident Plane Wave"; pp2 = "2Layers Thicknesses"; pp3 = "3Scatterer Properties"; @@ -27,7 +27,7 @@ DefineConstant[ thick_L_6 = {50 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]}, xsideg = {0 , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1}, - tag_geom = { 2 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar",8="cylindrical pillar"}}, + tag_geom = { 8 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar",8="cylindrical pillar"}}, rx = {60 , Name StrCat[pp3,"/1rx"]}, ry = {60 , Name StrCat[pp3,"/2ry"]}, rz = {130 , Name StrCat[pp3,"/3rz"]}, @@ -60,10 +60,10 @@ DefineConstant[ lc_scat = {lambda0/(5*paramaille) , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]}, PML_top = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]}, PML_bot = {lambda0, Name StrCat[pp5,"/5PML bot thickness [nm]"]}, - Nmax = {2 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]}, + Nmax = {1 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]}, refine_mesh_L_1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]}, - refine_mesh_L_2= {2 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]}, - refine_mesh_L_3= {6 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]}, + refine_mesh_L_2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]}, + refine_mesh_L_3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]}, refine_mesh_L_4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]}, refine_mesh_L_5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]}, refine_mesh_L_6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]}, diff --git a/DiffractionGratings/grating3D_parallel_Mmatrix.sh b/DiffractionGratings/grating3D_parallel_Mmatrix.sh index 25044f9..57b1005 100644 --- a/DiffractionGratings/grating3D_parallel_Mmatrix.sh +++ b/DiffractionGratings/grating3D_parallel_Mmatrix.sh @@ -2,50 +2,61 @@ export OPENBLAS_NUM_THREADS=1 export OMP_NUM_THREADS=1 export NPROC=96 -export nb_lam=50 -export nb_phi=50 -# export nb_lam=10 -# export nb_phi=10 + +export flag_angle_study="theta" +export loop_angle_max=50 +export fixed_angle=90 +export GRATING_CASE="pillar" + +# export flag_angle_study="phi" +# export loop_angle_max=360 +# export fixed_angle=50 +# export GRATING_CASE="half_ellipsoid" + +export nb_angle=2 +export nb_lam=2 export lambda_min=400 export lambda_max=1200 export FLAG_TOTAL=0 -export GRATING_CASE="half_ellipsoid" -export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL +export myDir="res_Matrix_nb_lam"$nb_lam"_nb_"$flag_angle_study$nb_angle"_total"$FLAG_TOTAL rm -r $myDir mkdir $myDir gmsh grating3D.geo -setstring test_case $GRATING_CASE -3 -o grating3D.msh myfunc() { - local mysubDir=$myDir/run_lam$1_phi$2_psi$3 + local mysubDir=$myDir/run_lam$1_$flag_angle_study$2_psi$3 mkdir $mysubDir cp grating3D.msh grating3D.pro grating3D_data_$GRATING_CASE.geo grating3D_data.geo grating3D_materials.pro $mysubDir - if (( $nb_lam == 1 )) - then - local lam=$(echo "scale=10;$lambda_min" | bc ) + local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc ) + if [ $flag_angle_study = "phi" ]; then + local phi=$(echo "scale=10;$loop_angle_max/($nb_angle-1)*$2" | bc ) + local theta=$(echo "scale=10;$fixed_angle" | bc ) else - local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc ) + local theta=$(echo "scale=10;$loop_angle_max/($nb_angle-1)*$2" | bc ) + local phi=$(echo "scale=10;$fixed_angle" | bc ) fi - local phi=$(echo "scale=10;360/($nb_phi)*$2" | bc ) local psi=$(echo "scale=10;90*$3" | bc ) cd $mysubDir - getdp grating3D.pro -pre helmholtz_vector -msh grating3D.msh -cal -pos postop_helmholtz_vector -petsc_prealloc 200 -setstring test_case $GRATING_CASE -setnumber lambda0 $lam -setnumber thetadeg 50 -setnumber phideg $phi -setnumber psideg $psi -setnumber FLAG_TOTAL $FLAG_TOTAL - if [ $3 -eq 0 ]; then cp res3D/rs.txt ../r_pin_sout_lam$1_phi$2.out ; fi - if [ $3 -eq 0 ]; then cp res3D/rp.txt ../r_pin_pout_lam$1_phi$2.out ; fi - if [ $3 -eq 1 ]; then cp res3D/rs.txt ../r_sin_sout_lam$1_phi$2.out ; fi - if [ $3 -eq 1 ]; then cp res3D/rp.txt ../r_sin_pout_lam$1_phi$2.out ; fi - if [ $3 -eq 0 ]; then cp res3D/ts.txt ../t_pin_sout_lam$1_phi$2.out ; fi - if [ $3 -eq 0 ]; then cp res3D/tp.txt ../t_pin_pout_lam$1_phi$2.out ; fi - if [ $3 -eq 1 ]; then cp res3D/ts.txt ../t_sin_sout_lam$1_phi$2.out ; fi - if [ $3 -eq 1 ]; then cp res3D/tp.txt ../t_sin_pout_lam$1_phi$2.out ; fi - cp res3D/eff_t1.txt ../eff_t1_lam$1_phi$2_psi$3.out - cp res3D/eff_r1.txt ../eff_r1_lam$1_phi$2_psi$3.out - cp res3D/eff_t2.txt ../eff_t2_lam$1_phi$2_psi$3.out - cp res3D/eff_r2.txt ../eff_r2_lam$1_phi$2_psi$3.out - cp res3D/temp-Q_scat.txt ../Q_scat_lam$1_phi$2_psi$3.out + getdp grating3D.pro -pre helmholtz_vector -msh grating3D.msh -cal -pos postop_helmholtz_vector -petsc_prealloc 200 -setstring test_case $GRATING_CASE -setnumber lambda0 $lam -setnumber thetadeg $theta -setnumber phideg $phi -setnumber psideg $psi -setnumber FLAG_TOTAL $FLAG_TOTAL + if [ $3 -eq 0 ]; then cp res3D/rs.txt ../r_pin_sout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 0 ]; then cp res3D/rp.txt ../r_pin_pout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 1 ]; then cp res3D/rs.txt ../r_sin_sout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 1 ]; then cp res3D/rp.txt ../r_sin_pout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 0 ]; then cp res3D/ts.txt ../t_pin_sout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 0 ]; then cp res3D/tp.txt ../t_pin_pout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 1 ]; then cp res3D/ts.txt ../t_sin_sout_lam$1_$flag_angle_study$2.out ; fi + if [ $3 -eq 1 ]; then cp res3D/tp.txt ../t_sin_pout_lam$1_$flag_angle_study$2.out ; fi + cp res3D/eff_t1.txt ../eff_t1_lam$1_$flag_angle_study$2_psi$3.out + cp res3D/eff_r1.txt ../eff_r1_lam$1_$flag_angle_study$2_psi$3.out + cp res3D/eff_t2.txt ../eff_t2_lam$1_$flag_angle_study$2_psi$3.out + cp res3D/eff_r2.txt ../eff_r2_lam$1_$flag_angle_study$2_psi$3.out + cp res3D/temp-Q_scat.txt ../Q_scat_lam$1_$flag_angle_study$2_psi$3.out cd ../.. rm -r $mysubDir } export -f myfunc -parallel -j $NPROC myfunc ::: $(seq 0 $(($nb_lam-1))) ::: $(seq 0 $(($nb_phi-1))) ::: 0 1 +parallel -j $NPROC myfunc ::: $(seq 0 $(($nb_lam-1))) ::: $(seq 0 $(($nb_angle-1))) ::: 0 1 + +python grating3D_postplot_Mmatrix.py -flag_angle_study $flag_angle_study -nb_angle $nb_angle -loop_angle_max $loop_angle_max -fixed_angle $fixed_angle -nb_lam $nb_lam -lambda_min $lambda_min -lambda_max $lambda_max diff --git a/DiffractionGratings/grating3D_postplot_Mmatrix.py b/DiffractionGratings/grating3D_postplot_Mmatrix.py index 5b7b783..c0b4590 100644 --- a/DiffractionGratings/grating3D_postplot_Mmatrix.py +++ b/DiffractionGratings/grating3D_postplot_Mmatrix.py @@ -2,6 +2,29 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import scipy.constants as scc +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("-flag_angle_study" , help="flag_angle_study" , type=str , default='phi') +parser.add_argument("-nb_angle" , help="nb_angle" , type=int , default=2) +parser.add_argument("-loop_angle_max" , help="loop_angle_max" , type=float , default=360) +parser.add_argument("-fixed_angle" , help="fixed_angle" , type=float , default=50) +parser.add_argument("-nb_lam" , help="nb_lam" , type=int , default=2) +parser.add_argument("-lambda_min" , help="lambda_min" , type=float , default=400) +parser.add_argument("-lambda_max" , help="lambda_max" , type=float , default=1200) +args = parser.parse_args() + +flag_angle_study = args.flag_angle_study +nb_angle = args.nb_angle +loop_angle_max = args.loop_angle_max +fixed_angle = args.fixed_angle +nb_lam = args.nb_lam +lambda_min = args.lambda_min +lambda_max = args.lambda_max + +FLAG_TOT = 0 + + plt.rcParams.update({"text.usetex": True, "font.family": "serif"}) def add_colorbar(mappable): @@ -23,14 +46,7 @@ def load_getdp_integral(fname): else: return temp[:,1]+1j*temp[:,2] -nb_lam = 50 -nb_angle = 51 -FLAG_TOT = 0 - -flag_study = 'theta' -# flag_study = 'phi' - -respath = 'res_Matrix_nb_lam%g_nb_%s%g_total0/'%(nb_lam,flag_study,nb_angle) +respath = 'res_Matrix_nb_lam%g_nb_%s%g_total0/'%(nb_lam,flag_angle_study,nb_angle) rpp = np.zeros((nb_lam,nb_angle),dtype=complex) rps = np.zeros((nb_lam,nb_angle),dtype=complex) @@ -47,13 +63,8 @@ effr1_sin = np.zeros((9,nb_lam,nb_angle),dtype=complex) effr2_sin = np.zeros((9,nb_lam,nb_angle),dtype=complex) Qscat_sin = np.zeros((nb_lam,nb_angle)) - -lambda_min = 400 -lambda_max = 1200 -angledeg_max = 50 - tab_lam = np.linspace(lambda_min,lambda_max,nb_lam) -tab_angle = np.linspace(0,angledeg_max,nb_angle) +tab_angle = np.linspace(0,loop_angle_max,nb_angle) tab_hnu = 2*scc.pi*scc.c/(tab_lam*1e-9/scc.hbar*scc.eV) M = np.zeros((4,4,nb_lam,nb_angle),dtype=complex) @@ -71,26 +82,20 @@ str_reftrans = 'r' for i in range(nb_lam): print(i) for j in range(nb_angle): - rpp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_pout_lam%g_%s%g.out'%(i,flag_study,j)) - rss[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_sout_lam%g_%s%g.out'%(i,flag_study,j)) - rps[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_pout_lam%g_%s%g.out'%(i,flag_study,j)) - rsp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_sout_lam%g_%s%g.out'%(i,flag_study,j)) - efft1_pin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi0.out'%(i,flag_study,j)) - efft2_pin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi0.out'%(i,flag_study,j)) - effr1_pin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi0.out'%(i,flag_study,j)) - effr2_pin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi0.out'%(i,flag_study,j)) - Qscat_pin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi0.out'%(i,flag_study,j))) - efft1_sin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi1.out'%(i,flag_study,j)) - efft2_sin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi1.out'%(i,flag_study,j)) - effr1_sin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi1.out'%(i,flag_study,j)) - effr2_sin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi1.out'%(i,flag_study,j)) - Qscat_sin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi1.out'%(i,flag_study,j))) - -# ## replicate data for phi=2\pi from phi=0 -# rpp[:,-1] = rpp[:,0] -# rss[:,-1] = rss[:,0] -# rps[:,-1] = rps[:,0] -# rsp[:,-1] = rsp[:,0] + rpp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_pout_lam%g_%s%g.out'%(i,flag_angle_study,j)) + rss[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_sout_lam%g_%s%g.out'%(i,flag_angle_study,j)) + rps[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_pout_lam%g_%s%g.out'%(i,flag_angle_study,j)) + rsp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_sout_lam%g_%s%g.out'%(i,flag_angle_study,j)) + efft1_pin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j)) + efft2_pin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j)) + effr1_pin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j)) + effr2_pin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j)) + Qscat_pin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j))) + efft1_sin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j)) + efft2_sin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j)) + effr1_sin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j)) + effr2_sin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j)) + Qscat_sin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j))) T1_sin = np.sum(efft1_sin,axis=0) T2_sin = np.sum(efft2_sin,axis=0) @@ -167,7 +172,7 @@ plt.savefig('BALANCE_new_code.jpg') fig, axes = plt.subplots(4, 4, subplot_kw=dict(projection='polar') ,figsize=(12,12)) flag_lam = False rlabel = r"$\hbar \nu$" if not flag_lam else r"$\lambda_0$" -anglelabel = r"$\varphi_0$" if flag_study=='phi' else r"$\theta_0$" +anglelabel = r"$\varphi_0$" if flag_angle_study=='phi' else r"$\theta_0$" which_r = tab_hnu if not flag_lam else tab_lam which_orig = 0 #if not flag_lam else 350 for i in range(4): @@ -198,7 +203,7 @@ for i in range(4): sm.set_array([]) ax.set_rorigin(which_orig) plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,wspace=0.35) -plt.savefig('Mmatrix_%s_%s.jpg'%(flag_study,str_reftrans)) +plt.savefig('Mmatrix_%s_%s.jpg'%(flag_angle_study,str_reftrans)) # plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0) -plt.show() +# plt.show() -- GitLab