I managed to write a SLEPc version that using EPSSolve. But the scale of spatial region is very tiny (about 1e-10), if I directly use EPSSolve, there will be error zero pivot. It looks like that SLEPc just ignore the problem.I have to rewrite th equation to scale.
How ever the original problem can be solved by ARPACK directly. Is there any way to scale the matrix in SLEPc automatically without rewrite the equation? I tried to search in SLEPc documentation but didn’t find any clue.
The code to solve original problem using SLEPc is( there will be error):
//Run: freefem++-mpi Schdroger3D.edp
real me=9.10938356e-31; //electron mass
real planckh=6.62607004e-34;
real hbar=1.0545718e-34;
real kgm2Ds2ToeV=6.24150636309e18;
real a0=5.2917721039e-11;//bohr radius
real qe=1.6021766208e-19;//electron charge
real varepsilon0=8.854187817e-12;
load "UMFPACK64"
load "msh3"
load "TetGen"
load "medit"
include "MeshSurface.idp"
load "PETSc"
macro dimension()3 // EOM
include "macro_ddm.idp"
//load "MUMPS_mpi"
int k=17;
int n=20;
real hs1 = a0; //mesh size on sphere
real hs2 = a0; //mesh size on sphere
func Pk=P1;
//正方形
int[int] nn = [n, n, n];
//Box: [-k*a0, k*a0]^3
real [int,int] Box = [[-k*a0, k*a0], [-k*a0, k*a0], [-k*a0, k*a0]]; //Bounding box: [xmin, xmax]
int [int,int] L = [[1, 2], [3, 4], [5, 6]]; //label: left,right, front, back, down, right
meshS ThH = SurfaceHex(nn, Box, L, 1);
meshS ThS1 = Sphere(1.5*a0, hs1, 7, 1);
meshS ThS2 = Sphere(5*a0, hs2, 8, 1);
meshS ThHS = ThH + ThS1 +ThS2;
//medit("Hex-Sphere", ThHS);
//Mesh
real[int] domain = [0, 0, 0, 1, hs1^3/6,
2.5*a0, 0,0, 2, hs2^3/6,
10*a0,0,0, 3, (2*k*a0/nn[1])^3/6];
mesh3 Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=3, regionlist=domain);
//Th=change(Th,fregion=nuTriangle%mpisize);
//medit("Cube with ball", Th);
fespace Vh(Th,Pk);
Vh Psi, v;
//real tilE=-2e-18;
varf Schrodinger(Psi,v) =
int3d(Th)( hbar^2/(2*me)*(dx(Psi)*dx(v) + dy(Psi)*dy(v) + dz(Psi)*dz(v)))
- int3d(Th)( qe^2/(4*pi*varepsilon0*sqrt(x^2+y^2+z^2))*Psi*v);
//+ on(1,2,3,4,5,6,Psi=0);
//- int3d(Th, mpirank)(tilE*Psi*v)
//+ on(1,2,3,4,5,6,Psi=0);
varf RHS(Psi,v) = int3d(Th)(Psi*v);
matrix locA = Schrodinger(Vh,Vh,solver="SPARSESOLVER");
matrix locB = RHS(Vh,Vh,solver="SPARSESOLVER",eps=1e-20);
Mat A;
createMat(Th, A, Pk)
A=locA;
Mat B(A, locB);
// 求解接近于sigma的特征值Ax=lambda*B*x，左边为A-sigma*B，右边为B。
int nev=10;
Vh[int] EigenPsi(nev);
real[int] EigenVal(nev);
int num = EPSSolve(A,B, sparams = " -st_type sinvert -eps_nev "+nev
+ " -eps_gen_hermitian"
+ " -eps_target -16.0",
values = EigenVal, vectors = EigenPsi);
//打印能级
for(int i = 0; i < nev; i++){
cout << "Eigenvalue["+i+"] = "<<EigenVal[i]
<<" , Energy level: " << EigenVal[i]*kgm2Ds2ToeV << " eV" << endl;
}
Scaled version is:
//verbosity=0;
include "constants.edp";
load "UMFPACK64"
load "msh3"
load "TetGen"
load "medit"
include "MeshSurface.idp";
load "PETSc"
macro dimension()3 // EOM
include "macro_ddm.idp";
//网格
real R=20; //Original range: [-20a0, 20a0]^3
real a=8.0; //2/(a0*k)=a
real k=2.0/(a0*a); //Scaled range : [-k*20*a0, k*20*a0]^3
real l=k*R*a0;
int n=25;
real R1=0.1*l;
real R2=0.4*l;
real hs1= 0.04*R1;
real hs2=0.1*R2;
func Pk=P1;
int[int] nn = [n, n, n];
//Box: [-k*a0, k*a0]^3
real [int,int] Box = [[-l, l], [-l, l], [-l, l]]; //Bounding box: [xmin, xmax]
int [int,int] L = [[1, 2], [3, 4], [5, 6]]; //label: left,right, front, back, down, right
meshS ThH = SurfaceHex(nn, Box, L, 1);
meshS ThS1 = Sphere(R1, hs1, 7, 1);
meshS ThS2 = Sphere(R2, hs2, 8, 1);
meshS ThHS = ThH + ThS1 +ThS2;
//medit("Hex-Sphere", ThHS);
//Mesh
real[int] domain = [0, 0, 0, 1, hs1^3/6,
(R1+R2)/2, 0,0, 2, hs2^3/6,
(R2+l)/2,0,0, 3, (2*l/nn[1])^3/6];
mesh3 Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=3, regionlist=domain);
medit("All", Th);
fespace Vh(Th,Pk);
Vh Psi, v;
varf Schrodinger(Psi,v) =
int3d(Th)(dx(Psi)*dx(v) + dy(Psi)*dy(v) + dz(Psi)*dz(v) - a/sqrt(x^2+y^2+z^2)*Psi*v )
//- int3d(Th)(tilE*Psi*v)
//+ on(1,2,3,4,5,6,Psi=0)
;
varf RHS(Psi,v) = int3d(Th)(Psi*v);
matrix locA = Schrodinger(Vh,Vh,solver="SPARSESOLVER");
matrix locB = RHS(Vh,Vh,solver="SPARSESOLVER",eps=1e-20);
Mat A;
createMat(Th, A, Pk)
A=locA;
Mat B(A, locB);
// 求解接近于sigma的特征值Ax=lambda*B*x，左边为A-sigma*B，右边为B。
int nev=10;
Vh[int] EigenPsi(nev);
real[int] EigenVal(nev);
int num = EPSSolve(A,B, sparams = " -st_type sinvert -eps_nev "+nev
+ " -eps_gen_hermitian"
+ " -eps_target -16.0",
values = EigenVal, vectors = EigenPsi);
for(int i = 0; i < nev; i++){
cout << "Eigenvalue["+i+"] = "<<EigenVal[i]
<<" , Energy level: " << EigenVal[i]*hbar^2*k^2/(2*me)*kgm2Ds2ToeV << " eV" << endl;
}