From 708bc2431e5b761e51b47ca6f0da663cc5b01021 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Tue, 27 Sep 2022 07:36:00 +0200
Subject: [PATCH] Bermudez PML

---
 DiffractionGratings/grating3D.pro             | 40 ++++++++++++-------
 .../grating3D_data_2Dlamellar.geo             |  1 +
 DiffractionGratings/grating3D_data_bisin.geo  |  3 +-
 .../grating3D_data_checker.geo                |  3 +-
 DiffractionGratings/grating3D_data_conv.geo   |  3 +-
 .../grating3D_data_halfellipsoid.geo          |  4 +-
 DiffractionGratings/grating3D_data_hole.geo   |  3 +-
 .../grating3D_data_pyramid.geo                |  1 +
 DiffractionGratings/grating3D_data_skew.geo   |  3 +-
 .../grating3D_data_solarcell.geo              |  2 +-
 DiffractionGratings/grating3D_data_torus.geo  |  3 +-
 11 files changed, 42 insertions(+), 24 deletions(-)

diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index b9396c3..8bebf20 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -125,27 +125,37 @@ Function{
 	a_pml           = 1.;
 	b_pml           = 1.;
   // bermu
-  n1[]    = Sqrt[epsr1[]];
-  n2[]    = Sqrt[epsr2[]];
-  k1norm[]= k0*n1[];
-  k2norm[]= k0*n2[];
-  Zmax = PML_top_hh;
-  Zmin = hh_L_6;
-  Damp_pml_top[] = 1/Abs[Zmax + PML_top - Z[]];
-  Damp_pml_bot[] = 1/Abs[Zmin - PML_bot - Z[]];
-  // DampingProfileX[] = 1/(Xmax + SizePMLX - Fabs[X[]]) - 1/(SizePMLX);
-
-  sx            = 1.;
-	sz[]          = Complex[a_pml,b_pml];
-	sz_bermutop[] = Complex[1,Damp_pml_top[]/k1norm[]];
-	sz_bermubot[] = Complex[1,Damp_pml_bot[]/k2norm[]];
-	sy            = 1.;
+  n1[]     = Sqrt[epsr1[]];
+  n2[]     = Sqrt[epsr2[]];
+  k1norm[] = k0*n1[];
+  k2norm[] = k0*n2[];
+
+  Zmax     = PML_top_hh;
+  Zmin     = hh_L_6;
+  Damp_pml_top[] = 1/(Zmax + PML_top - Fabs[Z[]]) - 1/(PML_top);
+  Damp_pml_bot[] = 1/(Zmin + PML_top - Fabs[Z[]]) - 1/(PML_bot);
+  Sigma_top[] = 0.5*(Damp_pml_top[] + Fabs[Damp_pml_top[]]);
+  Sigma_bot[] = 0.5*(Damp_pml_bot[] + Fabs[Damp_pml_bot[]]);
+
+	If (PML_TYPE==0)
+    sz[]          = Complex[a_pml,b_pml];
+	Else
+    sz[PMLtop] = Complex[1,Damp_pml_top[]/k1norm[]];
+	  sz[PMLbot] = Complex[1,Damp_pml_bot[]/k2norm[]];
+  EndIf
+  sx = 1.;
+	sy = 1.;
 
 	epsr[PMLtop]  = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 	mur[PMLtop]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
   epsr[PMLbot]  = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
   mur[PMLbot]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 
+  // epsr[PMLtop]  = Re[epsr1[]]*TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
+	// mur[PMLtop]   =             TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
+  // epsr[PMLbot]  = Re[epsr2[]]*TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+  // mur[PMLbot]   =             TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+
   // epsr[PMLbot] = epsr2[];
   // mur[PMLbot]  = TensorDiag[1,1,1];
 
diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
index 27dc260..ed3c4a7 100644
--- a/DiffractionGratings/grating3D_data_2Dlamellar.geo
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.geo
@@ -62,6 +62,7 @@ DefineConstant[
     refine_mesh_L_5= {0.3       , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {0.3       , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
     
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 0    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_bisin.geo b/DiffractionGratings/grating3D_data_bisin.geo
index 65840b2..216ef0b 100644
--- a/DiffractionGratings/grating3D_data_bisin.geo
+++ b/DiffractionGratings/grating3D_data_bisin.geo
@@ -62,7 +62,8 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
     FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
index 169cb1f..2e43fc2 100644
--- a/DiffractionGratings/grating3D_data_checker.geo
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -62,7 +62,8 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
     FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_conv.geo b/DiffractionGratings/grating3D_data_conv.geo
index 1e227c4..ebe61a7 100644
--- a/DiffractionGratings/grating3D_data_conv.geo
+++ b/DiffractionGratings/grating3D_data_conv.geo
@@ -63,7 +63,8 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets  = {0       , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-        
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = { 10 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1  , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
     FlagOutEtotCuts    = { 1  , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
index 04cd83c..e27f2d1 100644
--- a/DiffractionGratings/grating3D_data_halfellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -7,8 +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       = {400  , Min 400, Max 800, Step 200, Name StrCat[pp1, "0wavelength [nm]"] , Loop 1},
+    lambda0       = {495  , 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]"]},
@@ -63,6 +62,7 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
         
     InterpSampling     = { 10   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 0    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
index d057124..32bf712 100644
--- a/DiffractionGratings/grating3D_data_hole.geo
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -68,7 +68,8 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
     FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
index fb97f52..579361b 100644
--- a/DiffractionGratings/grating3D_data_pyramid.geo
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -62,6 +62,7 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
 
     InterpSampling     = { 50   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_skew.geo b/DiffractionGratings/grating3D_data_skew.geo
index 7ceee28..2bf61c7 100644
--- a/DiffractionGratings/grating3D_data_skew.geo
+++ b/DiffractionGratings/grating3D_data_skew.geo
@@ -63,7 +63,8 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-        
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = { 10   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
     FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_solarcell.geo b/DiffractionGratings/grating3D_data_solarcell.geo
index 7c26bb7..ccfd3c9 100644
--- a/DiffractionGratings/grating3D_data_solarcell.geo
+++ b/DiffractionGratings/grating3D_data_solarcell.geo
@@ -63,7 +63,7 @@ DefineConstant[
     refine_mesh_L_5= {2      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {2      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
     
     InterpSampling     = { 20   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
index 6435345..bc27227 100644
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ b/DiffractionGratings/grating3D_data_torus.geo
@@ -62,7 +62,8 @@ DefineConstant[
     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 [-]"]},
     FlagLinkFacets = {1      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = {15   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
     FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
-- 
GitLab