diff --git a/NonLinearEVP/NonLinearEVP.pro b/NonLinearEVP/NonLinearEVP.pro
new file mode 100644
index 0000000000000000000000000000000000000000..8c8ea73daaf4dd4ad5b83ea6e4399ecc084d2e08
--- /dev/null
+++ b/NonLinearEVP/NonLinearEVP.pro
@@ -0,0 +1,540 @@
+///////////////////////////////////
+//// Author : Guillaume Demesy ////
+//// .pro ////
+///////////////////////////////////
+
+Include "NonLinearEVP_data.geo";
+
+Group {
+ // Domains
+ Om_1 = Region[100] ; //in (dispersive)
+ Om_2inter = Region[101]; // out (non dispersive)
+ pmltop = Region[102]; // pml top
+ pmlbot = Region[103]; // pml top
+ Om = Region[{Om_1,Om_2inter,pmltop,pmlbot}];
+ Om_nopml = Region[{Om_1,Om_2inter}];
+ Om_pml = Region[{pmltop,pmlbot}];
+ Om_2 = Region[{Om_pml,Om_2inter}]; // non dispersive
+
+ // Boundaries
+ SurfBlochLeft = Region[1001];
+ SurfBlochRight = Region[1002];
+ Surfpml = Region[{1003,1004}];
+ Bound = Region[{1005}];
+ PrintPoint = Region[10000];
+}
+
+Function{
+ I[] = Complex[0,1];
+
+ siwt = -1.;
+ a_pml = 1.;
+ b_pml = -siwt;
+
+ sx[Om_pml] = 1.;
+ sy[Om_pml] = Complex[a_pml,b_pml];
+ sx[Om_nopml] = 1.;
+ sy[Om_nopml] = 1.;
+
+ sz = 1;
+
+ PML_Tensor[] = TensorDiag[sz*sy[]/sx[],sx[]*sz/sy[],sx[]*sy[]/sz];
+ mur[] = TensorDiag[1,1,1]*PML_Tensor[];
+ epsr_nod[] = TensorDiag[1,1,1]*PML_Tensor[];
+
+ term_o3_3cc_1[] = Complex[-eps_oo_1,0];
+ term_o3_2cc_1[] = Complex[0,-eps_oo_1*gam_1];
+ term_o3_1cc_1[] = Complex[om_d_1^2,0];
+ term_o3_1rr_1[] = cel^2*Complex[1,0];
+ term_o3_0rr_1[] = cel^2*Complex[0,gam_1];
+ term_o3_3cc_2[] = -1*epsr_nod[];
+ term_o3_2cc_2[] = -1*epsr_nod[]*Complex[0,gam_2];
+ term_o3_1cc_2[] = Complex[om_d_2^2,0];
+ term_o3_1rr_2[] = cel^2*Complex[1,0];
+ term_o3_0rr_2[] = cel^2*Complex[0,gam_2];
+ term_o3_3cc[Om_1] = term_o3_3cc_1[];
+ term_o3_2cc[Om_1] = term_o3_2cc_1[];
+ term_o3_1cc[Om_1] = term_o3_1cc_1[];
+ term_o3_1rr[Om_1] = term_o3_1rr_1[];
+ term_o3_0rr[Om_1] = term_o3_0rr_1[];
+ term_o3_3cc[Om_2] = term_o3_3cc_2[];
+ term_o3_2cc[Om_2] = term_o3_2cc_2[];
+ term_o3_1cc[Om_2] = term_o3_1cc_2[];
+ term_o3_1rr[Om_2] = term_o3_1rr_2[];
+ term_o3_0rr[Om_2] = term_o3_0rr_2[];
+
+
+ dephx[] = Complex[ Cos[kx*a_lat] , Sin[kx*a_lat]];
+ // dephy[] = Complex[ Cos[ky*a_lat] , Sin[ky*a_lat]];
+
+ // bidon
+ EigFilter[] = (Norm[$EigenvalueReal] > 1e-150);
+
+ //Auxialiary Field resolution coefficients
+ eps_inf[] = TensorDiag[1.,1.,1.]*eps_oo_1;
+ Om_D[] = TensorDiag[1,1,1]*Sqrt[om_d_1^2/( 2. )];
+ Gamma_d[] = TensorDiag[1,1,1]*gam_1;
+}
+
+Constraint {
+ {Name BlochX;
+ Case {
+ { Region SurfBlochRight; Type LinkCplx ; RegionRef SurfBlochLeft;
+ Coefficient dephx[]; Function Vector[$X-a_lat,$Y,$Z] ;
+ }
+ }
+ }
+
+ {Name Dirichlet; Type Assign;
+ Case {
+ { Region Surfpml ; Value 0.; }
+ }
+ }
+
+ // div condition
+ { Name Constraint_cInt; Type Assign;
+ Case {
+ }
+ }
+ { Name Constraint_Int; Type Assign;
+ Case {
+ { Region Bound; Value 0.; }
+ }
+ }
+}
+
+
+
+Jacobian {
+ {Name JVol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name JSur ;
+ Case{
+ { Region All ; Jacobian Sur ; }
+ }
+ }
+ { Name JLin ;
+ Case{
+ { Region All ; Jacobian Lin ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hgrad_perp; Type Form1P;
+ BasisFunction {
+ { Name un; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[All]; }
+ { Name un2; NameOfCoef un2; Function BF_PerpendicularEdge_2E; Support Region[Om]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name un3; NameOfCoef un3; Function BF_PerpendicularEdge_2F; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name un4; NameOfCoef un4; Function BF_PerpendicularEdge_3E; Support Region[Om]; Entity EdgesOf[All]; }
+ { Name un5; NameOfCoef un5; Function BF_PerpendicularEdge_3F; Support Region[Om]; Entity FacetsOf[All]; }
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ If(flag_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un4; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ EndIf
+ }
+ }
+
+ { Name Hcurl; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om]; Entity EdgesOf[All]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ If(flag_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
+ }
+
+ { Name Hcurl_aD; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om_1]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om_1]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om_1]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om_1]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om_1]; Entity EdgesOf[All]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ }
+
+ { Name Hcurl_1; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Om_1,Bound}]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Om_1,Bound}]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ }
+ { Name Hcurl_2; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Om_2,Bound}]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Om_2,Bound}]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ If(flag_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
+ }
+ { Name Hcurl_l; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Bound]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Bound]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ // { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Bound]; Entity FacetsOf[All]; }
+ // { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Bound]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Bound]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint {
+ If(flag_o2==1)
+
+ EndIf
+ }
+ }
+
+
+}
+
+Formulation {
+ { Name modal_Aux_E; Type FemEquation;
+ Quantity {
+ { Name e ; Type Local; NameOfSpace Hcurl ;}
+ { Name eaD; Type Local; NameOfSpace Hcurl_aD;}
+ }
+ Equation {
+ Galerkin{ [ 1.* cel^2/mur[] * Dof{Curl e}, {Curl e} ] ; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -1 *-epsr_nod[] * Dof{e} , {e} ] ; Order 2 ; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ [ 1 * 2.*om_d_1/Sqrt[2] * Dof{e} , {eaD} ] ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ [ -I[] * gam_1 * Dof{eaD} , {eaD} ] ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -I[] * Om_D[] * Dof{eaD } , {e} ] ; Order 1 ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -I[] *-Dof{eaD} , {eaD} ] ; Order 1 ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ { Name modal_helmholtz_Lag_E; Type FemEquation;
+ Quantity {
+ { Name u_1 ; Type Local ; NameOfSpace Hcurl_1 ;}
+ { Name u_2 ; Type Local ; NameOfSpace Hcurl_2 ;}
+ { Name lambda; Type Local ; NameOfSpace Hcurl_l ;}
+ }
+ Equation {
+ Galerkin { [ 1.* term_o3_0rr_1[]/mur[]* Dof{Curl u_1},{Curl u_1}]; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1cc_1[] * Dof{u_1} ,{u_1} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1rr_1[]/mur[]* Dof{Curl u_1},{Curl u_1}]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* term_o3_2cc_1[] * Dof{u_1} ,{u_1} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ I[] * term_o3_3cc_1[] * Dof{u_1} ,{u_1} ]; Order 3; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { [ 1. * term_o3_0rr_2[]/mur[]* Dof{Curl u_2},{Curl u_2}]; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1cc_2[] * Dof{u_2} ,{u_2} ]; Order 1;In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1rr_2[]/mur[]* Dof{Curl u_2},{Curl u_2}]; Order 1;In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1 * term_o3_2cc_2[] * Dof{u_2} ,{u_2} ]; Order 2;In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ I[] * term_o3_3cc_2[] * Dof{u_2} ,{u_2} ]; Order 3;In Om_2; Jacobian JSur; Integration Int_1;}
+
+ Galerkin { [ 1. * term_o3_0rr_1[] * 1e8 * Dof{lambda} , {u_1 } ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin { [ 1. * -term_o3_0rr_2[] * 1e8 * Dof{lambda} , {u_2 } ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin {Eig[ -I[] * term_o3_1rr_1[] * 1e8 * Dof{lambda} , {u_1 } ] ; Order 1; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin {Eig[ -I[] * -term_o3_1rr_2[] * 1e8 * Dof{lambda} , {u_2 } ] ; Order 1; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin { [ Dof{u_1} * 1e8 , {lambda} ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin { [ -Dof{u_2} * 1e8 , {lambda} ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ }
+ }
+ { Name modal_helmholtz_PEP_E; Type FemEquation;
+ Quantity {{ Name u ; Type Local; NameOfSpace Hcurl ;}}
+ Equation {
+ Galerkin { [ 1.* cel^2/mur[]*I[]*gam_1 * Dof{Curl u}, {Curl u}]; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* cel^2/mur[] * Dof{Curl u}, {Curl u}]; Order 1; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* om_d_1^2 * Dof{u} , {u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-1. * I[]*-eps_oo_1*gam_1 * Dof{u} , {u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[I[] * -eps_oo_1 * Dof{u} , {u} ]; Order 3; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-1. * -epsr_nod[]*I[]*gam_1 * Dof{u} , {u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[I[] * -epsr_nod[] * Dof{u} , {u} ]; Order 3; In Om_2; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name modal_helmholtz_PEP_E_nleig; Type FemEquation;
+ Quantity {{ Name u ; Type Local; NameOfSpace Hcurl ;}}
+ Equation {
+ Galerkin { Eig[ 1/mur[] * Dof{Curl u}, {Curl u} ]; Rational 1; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -epsr_nod[]/cel^2 * Dof{u}, {u} ]; Rational 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -epsr_nod[]/cel^2 * Dof{u}, {u} ]; Rational 3; In Om_2; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name modal_helmholtz_PEP_h; Type FemEquation;
+ Quantity {{ Name u ; Type Local; NameOfSpace Hgrad_perp ;}}
+ Equation {
+ Galerkin { [ 1.* 1/epsr_nod[] * -om_d_1^2 * Dof{Curl u}, {Curl u} ]; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* 1/epsr_nod[] * I[]*eps_oo_1*gam_1* Dof{Curl u}, {Curl u} ]; Order 1; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* 1/epsr_nod[] * eps_oo_1 * Dof{Curl u}, {Curl u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* 1/epsr_nod[] *I[]*gam_1 * Dof{Curl u}, {Curl u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* 1/epsr_nod[] * Dof{Curl u}, {Curl u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* mur[]/cel^2 * om_d_1^2 * Dof{u}, {u} ]; Order 2; In Om; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[I[] * -mur[]/cel^2 * I[]*eps_oo_1*gam_1* Dof{u}, {u} ]; Order 3; In Om; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ 1.* -mur[]/cel^2 * eps_oo_1 * Dof{u}, {u} ]; Order 4; In Om; Jacobian JSur; Integration Int_1;}
+ }
+ }
+}
+
+Resolution {
+ { Name res_PEP_E;
+ System{
+ { Name M1; NameOfFormulation modal_helmholtz_PEP_E ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1];EigenSolve[M1,neig,eig_target_re,eig_target_im];
+ // Print[M1];
+ }
+ }
+ { Name res_NEP_E;
+ System{{ Name M1; NameOfFormulation modal_helmholtz_PEP_E_nleig ; Type ComplexValue; }}
+ Operation{
+ GenerateSeparate[M1];
+ // // Fan rational (iw)
+ // (2*Pi*cel/lambda_vp)
+ EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[],
+ { {1}, {-eps_oo_1,gam_1*eps_oo_1, -om_d_1^2,0}, {-1,0,0} } ,
+ { {1}, {1,-gam_1}, {1} } ];
+ }
+ }
+
+ { Name res_Lag_E;
+ System{
+ { Name M1; NameOfFormulation modal_helmholtz_Lag_E ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im]; // Print[M1];
+ }
+ }
+ { Name res_PEP_h;
+ System{
+ { Name M1; NameOfFormulation modal_helmholtz_PEP_h ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1];
+ EigenSolve[M1,neig,eig_target_re,eig_target_im]; // Print[M1];
+ }
+ }
+
+ { Name res_Aux_E;
+ System{
+ { Name M1; NameOfFormulation modal_Aux_E ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[]];
+ }
+ }
+}
+
+PostProcessing {
+ { Name postpro_eigenvectors_PEP_E; NameOfFormulation modal_helmholtz_PEP_E;
+ Quantity {
+ { Name intnormu_pml ; Value { Integral { [ Norm[{u}] ] ; In Om_pml ; Jacobian JSur; Integration Int_1 ;} } }
+ { Name intnormu_nopml ; Value { Integral { [ Norm[{u}] ] ; In Om_nopml; Jacobian JSur; Integration Int_1 ;} } }
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_NEP_E; NameOfFormulation modal_helmholtz_PEP_E_nleig;
+ Quantity {
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_Lag_E; NameOfFormulation modal_helmholtz_Lag_E;
+ Quantity {
+ { Name u_1 ; Value { Local { [ {u_1} ] ; In Om_1; Jacobian JSur; } } }
+ { Name u_2 ; Value { Local { [ {u_2} ] ; In Om_2; Jacobian JSur; } } }
+ { Name lambda ; Value { Local { [ {lambda} ] ; In Bound ; Jacobian JLin; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_PEP_h; NameOfFormulation modal_helmholtz_PEP_h;
+ Quantity {
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_Aux_E; NameOfFormulation modal_Aux_E;
+ Quantity {
+ { Name e ; Value { Local { [ {e} ] ; In Om; Jacobian JSur; } } }
+ { Name eaD ; Value { Local { [ {eaD} ] ; In Om_1; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+}
+
+
+PostOperation {
+ { Name postop_PEP_E; NameOfPostProcessing postpro_eigenvectors_PEP_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u , OnElementsOf Om,File "u_PEP_E.pos", Name "Electric eigen-field (PEP-E)"];
+ EndIf
+ }
+ }
+ { Name postop_NEP_E; NameOfPostProcessing postpro_eigenvectors_NEP_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u , OnElementsOf Om,File "u_NEP_E.pos", Name "Electric eigen-field NEP-E"];
+ EndIf
+ }
+ }
+ { Name postop_Lag_E; NameOfPostProcessing postpro_eigenvectors_Lag_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u_1 , OnElementsOf Om_1 , File "u1_Lag_E.pos", Name "Electric eigen-field 1 (Lag-E)"];
+ Print [ u_2 , OnElementsOf Om_2 , File "u2_Lag_E.pos", Name "Electric eigen-field 2 (Lag-E)"];
+ Print [ lambda , OnElementsOf Bound, File "lam_Lag_E.pos", Name "Lagrange multiplier (Lag-E)"];
+ EndIf
+ }
+ }
+ { Name postop_PEP_h; NameOfPostProcessing postpro_eigenvectors_PEP_h ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u , OnElementsOf Om, File "u_PEP_h.pos", Name "Magnetic eigen-field (PEP-h)"];
+ EndIf
+ }
+ }
+ { Name postop_Aux_E; NameOfPostProcessing postpro_eigenvectors_Aux_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ eaD , OnElementsOf Om_1, File "u_Aux_EaD.pos", Name "Auxiliary Field Aux-E"];
+ Print [ e , OnElementsOf Om , File "u_Aux_E.pos", Name "E Aux-E"];
+ EndIf
+ }
+ }
+}
+
+// In the list of changes of PETSc 3.9
+// http://www.mcs.anl.gov/petsc/documentation/changes/39.html
+// MatSolverPackage is replaced with MatSolverType
+//
+// So you have to change in the command line option:
+// -st_pc_factor_mat_solver_package mumps
+// to:
+// -st_pc_factor_mat_solver_type mumps
+
+// fine tuning slecp options
+// -pep_ncv 600 (set size of the Krylov subspace)
+// -pep_mpd 600
+// -pep_toar_locking 0 ("do not take an EV for granted")
+
+
+eig_tol = 1e-6;
+slepc_options_nep_pol = " -nep_max_it 1000 -nep_type nleigs -nep_target_magnitude ";
+// Warning slepc > 3.8 : -nep_nleigs_rational => -nep_rational
+// slepc_options_nep_rat = Sprintf(" -nep_max_it 1000 -nep_type nleigs -nep_nleigs_rational -nep_target_magnitude -nep_tol %.10f ",eig_tol);
+slepc_options_nep_rat = Sprintf(" -nep_max_it 1000 -nep_type nleigs -nep_rational -nep_target_magnitude -nep_tol %.10f ",eig_tol);
+
+// Warning petsc > 3.9 st_pc_factor_mat_solver_package => st_pc_factor_mat_solver_type
+// slepc_options_st = "-solve -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_type mumps -petsc_prealloc 200 -pos";
+slepc_options_st = "-solve -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_package mumps -petsc_prealloc 200 -pos";
+
+slepc_options_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f ",eig_tol);
+
+If (flag_res==0)
+ which_res="Aux_E";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+If (flag_res==1)
+ which_res="PEP_E";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+If (flag_res==2)
+ which_res="NEP_E";
+ str_slepc_options=StrCat(slepc_options_nep_rat);
+EndIf
+
+If (flag_res==3)
+ which_res="Lag_E";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+If (flag_res==4)
+ which_res="PEP_h";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+// redirect converged eigenvalues at runtime to some textfile
+str_slepc_options = StrCat(str_slepc_options," -nep_monitor_all :temp_eigenvalues.txt ");
+
+DefineConstant[
+ R_ = {StrCat("res_",which_res), Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {Sprintf(StrCat("-solve -pos -v2 ",str_slepc_options,slepc_options_rg)), Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {StrCat("postop_",which_res), Name "GetDP/2PostOperationChoices", Visible 1}];