diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 0348ee7954691006089154c2aae08569146b65d8..b745ee53e591de6a91f0c9c9fd67886a1a72658e 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -22,15 +22,15 @@ myDir = "res3D/";
Group {
// SubDomains
- PMLbot = Region[1];
- L_6_temp = Region[2];
- L_5 = Region[3];
- L_4 = Region[4];
- L_3 = Region[5];
- L_2 = Region[6];
- L_1_temp = Region[7];
- PMLtop = Region[8];
- Scat = Region[9];
+ PMLbot = Region[1];
+ L_6 = Region[2];
+ L_5 = Region[3];
+ L_4 = Region[4];
+ L_3 = Region[5];
+ L_2 = Region[6];
+ L_1 = Region[7];
+ PMLtop = Region[8];
+ Scat = Region[9];
// Boundaries
SurfBlochXm = Region[101];
@@ -50,9 +50,6 @@ Group {
SurfExcludeFacets = Region[{SurfBloch}];
EndIf
- L_1 = Region[{L_1_temp,SurfIntTop}];
- L_6 = Region[{L_6_temp,SurfIntBot}];
-
Omega = Region[{PMLbot,L_6,L_5,L_4,L_3,L_2,L_1,PMLtop,Scat}];
Omega_nosource = Region[{PMLbot,L_6,L_1,PMLtop}];
Omega_source = Region[{Scat,L_2,L_3,L_4,L_5}];
@@ -209,14 +206,25 @@ Function{
Hr[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
Ht[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
+ E1[SurfIntTop] = Ei[]+Er[];
E1[Omega_super] = Ei[]+Er[];
E1[Omega_subs] = Et[];
+ E1[SurfIntBot] = Et[];
+
+ E1d[SurfIntTop] = Er[];
E1d[Omega_super] = Er[];
E1d[Omega_subs] = Et[];
+ E1d[SurfIntBot] = Et[];
+
+ H1[SurfIntTop] = Hi[]+Hr[];
H1[Omega_super] = Hi[]+Hr[];
H1[Omega_subs] = Ht[];
+ H1[SurfIntBot] = Ht[];
+
+ H1d[SurfIntTop] = Hr[];
H1d[Omega_super] = Hr[];
H1d[Omega_subs] = Ht[];
+ H1d[SurfIntBot] = Ht[];
source_vol_scat[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
source_surf_tot[] = -2*I[]*om0*mu0*Vector[0,0,1] /\ Hi[];
@@ -292,12 +300,12 @@ Integration {
FunctionSpace {
{ Name Hcurl; Type Form1;
BasisFunction {
- { Name sn; NameOfCoef un; Function BF_Edge; Support Region[{Omega}]; Entity EdgesOf[All]; }
- { Name sn2; NameOfCoef un2; Function BF_Edge_2E;Support Region[{Omega}]; Entity EdgesOf[All]; }
+ { Name sn; NameOfCoef un; Function BF_Edge; Support Region[{Omega,Gama}]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E;Support Region[{Omega,Gama}]; Entity EdgesOf[All]; }
If(oi==2)
- { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b; Support Region[Omega]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
- { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c; Support Region[Omega]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
- { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Omega]; Entity EdgesOf[Omega, Not SurfExcludeFacets]; }
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b; Support Region[{Omega,Gama}]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c; Support Region[{Omega,Gama}]; Entity FacetsOf[Omega, Not SurfExcludeFacets]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Omega,Gama}]; Entity EdgesOf[Omega, Not SurfExcludeFacets]; }
EndIf
}
Constraint {
diff --git a/DiffractionGratings/grating3D_data_half_ellipsoid.geo b/DiffractionGratings/grating3D_data_half_ellipsoid.geo
index 12f3e2a5d64aa8190d4f22fa045371561cb9e520..0f3f7881b0587561ffca9aa9c5eb9b020bb1cc81 100644
--- a/DiffractionGratings/grating3D_data_half_ellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_half_ellipsoid.geo
@@ -1,4 +1,4 @@
-nm = 1000;
+nm = 1;
pp1 = "1Incident Plane Wave";
pp2 = "2Layers Thicknesses";
pp3 = "3Scatterer Properties";
@@ -7,7 +7,7 @@ pp5 = "5Computational Parameters";
pp6 = "6Output";
DefineConstant[
FLAG_TOTAL = {0 , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
- lambda0 = {495 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ lambda0 = {1000 , Name StrCat[pp1,"/1lambda0 [nm]"]},
thetadeg = {40 , Name StrCat[pp1,"/2theta0 [deg]"]},
phideg = {36 , Name StrCat[pp1,"/3phi0 [deg]"]},
psideg = {72 , Name StrCat[pp1,"/4psi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_parallel_Mmatrix.sh b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
index 887d864165e8fa66fed6838ae31e15c1725e05ba..a0d966472363960e3cc3f55c4db9007d6469bdcd 100644
--- a/DiffractionGratings/grating3D_parallel_Mmatrix.sh
+++ b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
@@ -2,12 +2,10 @@
export OPENBLAS_NUM_THREADS=1
export OMP_NUM_THREADS=1
export NPROC=96
-# export nb_lam=1
-# export nb_phi=1
-export nb_lam=100
-export nb_phi=59
-export lambda_min=400
-export lambda_max=1600
+export nb_phi=72
+export nb_lam=50
+export lambda_min=495
+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
@@ -20,11 +18,11 @@ myfunc() {
local mysubDir=$myDir/run_lam$1_phi$2_psi$3
mkdir $mysubDir
cp grating3D.msh grating3D.pro grating3D_data_$GRATING_CASE.geo grating3D_data.geo grating3D_materials.pro $mysubDir
- local lam=$(echo "scale=10;400+400/($nb_lam-1)*$1" | bc )
+ local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc )
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
+ 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 55 -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
diff --git a/DiffractionGratings/grating3D_postplot_Mmatrix.py b/DiffractionGratings/grating3D_postplot_Mmatrix.py
index e224f275fb73d3b4ab5e6b8a693b4c4f5a0c1c98..e6b2b8dcad991b3a3700c9cd2f0f58e652e38fe9 100644
--- a/DiffractionGratings/grating3D_postplot_Mmatrix.py
+++ b/DiffractionGratings/grating3D_postplot_Mmatrix.py
@@ -23,10 +23,10 @@ def load_getdp_integral(fname):
else:
return temp[:,1]+1j*temp[:,2]
-nb_lam = 100
-nb_phi = 60
+nb_lam = 50
+nb_phi = 73
FLAG_TOT = 0
-respath = 'res_Matrix_nb_lam100_nb_phi59_total0/'
+respath = 'res_Matrix_nb_lam%g_nb_phi%g_total0/'%(nb_lam,nb_phi-1)
rpp = np.zeros((nb_lam,nb_phi),dtype=complex)
rps = np.zeros((nb_lam,nb_phi),dtype=complex)
rsp = np.zeros((nb_lam,nb_phi),dtype=complex)
@@ -43,7 +43,7 @@ effr2_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
Qscat_sin = np.zeros((nb_lam,nb_phi))
-tab_lam = np.linspace(350,1600,nb_lam)
+tab_lam = np.linspace(495,1200,nb_lam)
tab_phi = np.linspace(0,360,nb_phi)
tab_hnu = 2*scc.pi*scc.c/(tab_lam*1e-9/scc.hbar*scc.eV)
M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
@@ -91,6 +91,10 @@ T1_pin = np.sum(efft1_pin,axis=0)
T2_pin = np.sum(efft2_pin,axis=0)
R1_pin = np.sum(effr1_pin,axis=0)
R2_pin = np.sum(effr2_pin,axis=0)
+T00_sin = efft2_sin[4,:,:]
+R00_sin = effr2_sin[4,:,:]
+T00_pin = efft2_pin[4,:,:]
+R00_pin = effr2_pin[4,:,:]
TOT1_pin = T1_pin+R1_pin+Qscat_pin
TOT2_pin = T2_pin+R2_pin+Qscat_pin
TOT1_sin = T1_sin+R1_sin+Qscat_sin
@@ -139,6 +143,18 @@ for form in range(2):
axes[2,1].title.set_text('effr1 pin')
plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
+fig, axes = plt.subplots(3, 2, figsize=(12,12))
+L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
+ax=axes[0,0];zplot = ax.pcolormesh(L,P,T00_sin.real );add_colorbar(zplot);ax.title.set_text('T00 sin')
+ax=axes[1,0];zplot = ax.pcolormesh(L,P,R00_sin.real );add_colorbar(zplot);ax.title.set_text('R00 sin')
+ax=axes[2,0];zplot = ax.pcolormesh(L,P,Qscat_sin.real);add_colorbar(zplot);ax.title.set_text('Abs sin')
+ax=axes[0,1];zplot = ax.pcolormesh(L,P,T00_pin.real );add_colorbar(zplot);ax.title.set_text('T00 pin')
+ax=axes[1,1];zplot = ax.pcolormesh(L,P,R00_pin.real );add_colorbar(zplot);ax.title.set_text('R00 pin')
+ax=axes[2,1];zplot = ax.pcolormesh(L,P,Qscat_pin.real);add_colorbar(zplot);ax.title.set_text('Abs pin')
+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$"
@@ -153,9 +169,9 @@ for i in range(4):
ax.text(0. , 1. , r"$M_{%d%d}$"%(i+1,j+1) , fontsize=16 , transform=ax.transAxes)
ax.set_xticks([])
ax.set_yticks([])
- cbar = plt.colorbar(sm, ax=ax, fraction=0.046, pad=0.04)
- cbar.ax.locator_params(nbins=5)
- cbar.ax.tick_params(labelsize=14)
+ # cbar = plt.colorbar(sm, ax=ax, fraction=0.046, pad=0.04)
+ # cbar.ax.locator_params(nbins=5)
+ # cbar.ax.tick_params(labelsize=14)
else:
p00 = ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0]-1,cmap=plt.cm.bwr,vmin=-1,vmax=1)
ax.xaxis.label.set_color('C0') #setting up X-axis label color to yellow