From c910736050abd42628af4e7cec1ce52ad7ff4435 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Mon, 10 Sep 2018 17:39:53 +0200
Subject: [PATCH] (iw) is the eigenvalue => clearer if formulations do not
involve I[]
---
NonLinearEVP/NonLinearEVP.pro | 97 ++++++++++++-----------------------
1 file changed, 32 insertions(+), 65 deletions(-)
diff --git a/NonLinearEVP/NonLinearEVP.pro b/NonLinearEVP/NonLinearEVP.pro
index 137232a..de28e2b 100644
--- a/NonLinearEVP/NonLinearEVP.pro
+++ b/NonLinearEVP/NonLinearEVP.pro
@@ -25,7 +25,7 @@ Group {
}
Function{
- I[] = Complex[0,1];
+ I[] = Complex[0,1];
siwt = -1.;
a_pml = 1.;
@@ -65,15 +65,14 @@ Function{
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);
+ // dummy filter (better to use slepc to filter)
+ EigFilter[] = (Norm[$EigenvalueReal] > 1e-20);
//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;
+ Om_D[] = TensorDiag[1.,1.,1.]*om_d_1/Sqrt[2.];
+ Gamma_d[] = TensorDiag[1.,1.,1.]*gam_1;
}
Constraint {
@@ -90,7 +89,6 @@ Constraint {
{ Region Surfpml ; Value 0.; }
}
}
-
// div condition
{ Name Constraint_cInt; Type Assign;
Case {
@@ -103,8 +101,6 @@ Constraint {
}
}
-
-
Jacobian {
{Name JVol ;
Case {
@@ -207,7 +203,6 @@ FunctionSpace {
{ 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
}
}
@@ -219,7 +214,6 @@ FunctionSpace {
{ 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 {
@@ -242,8 +236,6 @@ FunctionSpace {
{ 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
}
@@ -264,12 +256,12 @@ Formulation {
{ 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;}
+ Galerkin{ [ cel^2/mur[] * Dof{Curl e}, {Curl e} ] ; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ epsr_nod[] * Dof{e} , {e} ] ; Order 2 ; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ [ Sqrt[2]*om_d_1 * Dof{e} , {eaD} ] ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ [ gam_1 * Dof{eaD} , {eaD} ] ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ Om_D[] * Dof{eaD } , {e} ] ; Order 1 ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -Dof{eaD} , {eaD} ] ; Order 1 ; In Om_1 ; Jacobian JSur; Integration Int_1;}
}
}
{ Name modal_helmholtz_Lag_E; Type FemEquation;
@@ -289,7 +281,7 @@ Formulation {
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;}
-
+ // The actual value of the lagrange multiplier is irrelevant => scaled to balance the eigenvectors
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;}
@@ -301,13 +293,13 @@ Formulation {
{ 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;}
+ Galerkin { [-cel^2/mur[]*gam_1*Dof{Curl u},{Curl u}]; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ cel^2/mur[] *Dof{Curl u},{Curl u}]; Order 1; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ om_d_1^2 *Dof{u} ,{u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-eps_oo_1*gam_1 *Dof{u} ,{u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ eps_oo_1 *Dof{u} ,{u} ]; Order 3; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-epsr_nod[]*gam_1 *Dof{u} ,{u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ epsr_nod[] *Dof{u} ,{u} ]; Order 3; In Om_2; Jacobian JSur; Integration Int_1;}
}
}
{Name modal_helmholtz_NEP_E; Type FemEquation;
@@ -321,14 +313,14 @@ Formulation {
{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;}
+ Galerkin { [-1/epsr_nod[]* om_d_1^2 * Dof{Curl u}, {Curl u} ]; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ 1/epsr_nod[]*eps_oo_1*gam_1* Dof{Curl u}, {Curl u} ]; Order 1; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-1/epsr_nod[]*eps_oo_1 * Dof{Curl u}, {Curl u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ 1/epsr_nod[]*gam_1 * Dof{Curl u}, {Curl u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-1/epsr_nod[] * Dof{Curl u}, {Curl u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-mur[]/cel^2 *om_d_1^2 * Dof{u}, {u} ]; Order 2; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ mur[]/cel^2 *eps_oo_1*gam_1* Dof{u}, {u} ]; Order 3; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-mur[]/cel^2 *eps_oo_1 * Dof{u}, {u} ]; Order 4; In Om ; Jacobian JSur; Integration Int_1;}
}
}
}
@@ -383,23 +375,6 @@ Resolution {
GenerateSeparate[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[]];
}
}
-
- // { Name all;
- // System{
- // { Name M1; NameOfFormulation modal_Aux_E ; Type ComplexValue; }
- // { Name M2; NameOfFormulation modal_helmholtz_PEP_E ; Type ComplexValue; }
- // { Name M3; NameOfFormulation modal_helmholtz_NEP_E ; Type ComplexValue; }
- // { Name M4; NameOfFormulation modal_helmholtz_Lag_E ; Type ComplexValue; }
- // { Name M5; NameOfFormulation modal_helmholtz_PEP_h ; Type ComplexValue; }
- // }
- // Operation{
- // GenerateSeparate[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[]];
- // GenerateSeparate[M2]; EigenSolve[M2,neig,eig_target_re,eig_target_im,EigFilter[]];
- // GenerateSeparate[M3]; EigenSolve[M3,neig,eig_target_re,eig_target_im,EigFilter[]];
- // GenerateSeparate[M4]; EigenSolve[M4,neig,eig_target_re,eig_target_im,EigFilter[]];
- // GenerateSeparate[M5]; EigenSolve[M5,neig,eig_target_re,eig_target_im,EigFilter[]];
- // }
- // }
}
PostProcessing {
@@ -498,15 +473,12 @@ PostOperation {
}
eig_tol = 1e-6;
-// slepc_options_nep_rat = Sprintf(" -nep_max_it 100 -nep_type nleigs -nep_rational -nep_tol %.10f -nep_view ",eig_tol);
-// slepc_options_nep_rat = Sprintf(" -nep_type nleigs -nep_rational -nep_view ");
slepc_options_nep_rat = Sprintf(" -nep_view ");
-//
-// Warning petsc > 3.9 st_pc_factor_mat_solver_package => st_pc_factor_mat_solver_type
-// TODO all this should now be default
-slepc_options_st = " -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_type mumps ";
-
-slepc_options_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f -pep_view ",eig_tol);
+// // These slepc options are now by default
+// slepc_options_st = " -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_type mumps ";
+// slepc_options_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f -pep_view ",eig_tol);
+slepc_options_st = " ";
+slepc_options_pep = Sprintf(" -pep_view ");
If (flag_res==0)
which_res="Aux_E";
@@ -533,11 +505,6 @@ If (flag_res==4)
str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
EndIf
-// If (flag_res==5)
-// which_res="all";
-// 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 ");
--
GitLab