I’m still trying to learn both FF and numerical methods details - initially I was going to do some
EM problems but wanted to get a fun 3D geometry and had meshing issues.
So, I backed up to 2D QM issue of interest . This code has been hacked up terribly
for learning but it found the lowest eigenvalue for sch eqn by complex time
stepping ( see the reference) which is similar apparently to the “power method”
to find largest eigenvalue but this finds the ground state. It seems to work well
by projecting out the prior eigenvectors. Curious if anyone cares to comment
on any utility of this approach or how it should be done and similar with the
FF coding. Thanks.
// run with MPI: ff-mpirun -np 4 script.edp
// NBPROC 4
//
// code borrowed from examples/ffddm/heat-torus-3d-surf.edp
/*
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ schtsmatR.edp -dt .001 -vc .01 -nx 200 -ny 200
@software{mjm_cpp_$PROJ,
author = {Michael J Marchywka},
title = {$PROJ},
abstract=(),
institution={},
license={Knowledge sir should be free to all },
publisher={Mike Marchywka},
email={marchywka@hotmail.com},
authorid={orcid.org/0000-0001-9237-455X},
filename = {$PROJ},
url = {},
version = {0.0.0},
date-started = {$DATE}
}
Complex time steps the schrodeinger eqn to get lower eigenvalue
and project out those already found.
[1]Solving the Schr\"{o}dinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients, Bader , Philipp[ ...] Casas , Fernando; The Journal of Chemical Physics. 2013 09/28/2013
https://arxiv.org/pdf/1304.6845.pdf
[2]How the inter-electronic potential Ans\"{a}tze affect the bound state solutions of a planar two-electron quantum dot model, Caruso , F.[ ...] Silveira , F.; Physica E: Low-dimensional Systems and Nanostructures. 2019 2019
https://arxiv.org/pdf/1809.08326.pdf
Did not do too bad but some out of order,
> x=1:25; k=1; c=pi*pi/100/100; for( i in 1:5) { for ( j in 1:5 ) { w=i*i+j*j; x[k]=c*w;k=k+1; }} ; sort(x)
[1] 0.001973921 0.004934802 0.004934802 0.007895684 0.009869604 0.009869604
[7] 0.012830486 0.012830486 0.016778327 0.016778327 0.017765288 0.019739209
[13] 0.019739209 0.024674011 0.024674011 0.025660971 0.025660971 0.028621853
[19] 0.028621853 0.031582734 0.033556655 0.033556655 0.040465378 0.040465378
[25] 0.049348022
>
2854 FreeFem++ schpro2.edp -dt .001 -vc 0 -nx 200 -ny 200
energy 0 0.00197335
energy 1 0.00493256
energy 2 0.00986318
energy 3 0.00493234
energy 4 0.00788999
energy 5 0.00985347
energy 6 0.0128064
energy 7 0.0177122
energy 8 0.0128016
energy 9 0.0167442
energy 10 0.0167405
energy 11 0.0196924
energy 12 0.025572
energy 13 0.0196778
marchywka@happy:/home/documents/cpp/proj/freefem/play$
*/
load "medit" // HPDDM plugin
include "macro_ddm.idp" // additional DDM functions
//https://doc.freefem.org/documentation/BEM.html
load "msh3"
// syntax errors due to forgotten macro are a huge problem
// and not clear from eror message !!!!!!!!!!!!!!!!!!!!
//macro init(i)i// EOM // scalar field initialization
macro Grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient
func Pk = P2; // finite element space
// box - to + or 2x in length
int szx = getARGV("-szx", 50); // refinement factor
int szy = getARGV("-szy", 50); // refinement factor
int nx =getARGV("-nx", 100);
int ny = getARGV("-ny", 100);
real dt = getARGV("-dt", 1e-4); // time step
int itermax = getARGV("-itermax", 10000); // number of iterations
int nstates = getARGV("-nstates", 20); // number of energies to keep
// error threshold to view and accept
real l2dun= getARGV("-l2dun", 1e-4); //
real l2view= getARGV("-l2view", 2e-4); //
// potential params
int sz2 = getARGV("-sz2", 10); // refinement factor
real vcenter= getARGV("-vc", 0.0); // V scale factor
cout<<vcenter<<endl;
real vx= getARGV("-vx", 1.0); // scale factor in x for V
real vy= getARGV("-vy", 1.0); // scale factoor in y for V
cout<<" szx "<<szx<<" szy "<<szy<<" sz2 "<<sz2<<endl;
// 2real R = 3, r=1, h=0.2;
//macro sum(v,xx) { real s=0; for(int i=0; i<xx[].n; ++i) s=s+xx[][i]; v= s; } // EOM
//func real a(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/sz; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1); if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm; }
func real ax(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/szx; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1); if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm; }
func real ay(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/szy; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1); if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm; }
//func fx= -szx + ax(x)*2*szx;
func fx= -szx + x*2*szx;
//func fx= -szx + 2*x*(.1*(x-.5)*(x-.5)+0*(y-.5)*(y-.5))*2*szx;
//func fx= -szx + 4*(x-.5)*((x-.5)*(x-.5)+0*(y-.5)*(y-.5))*2*szx*szx*szx;
//func fy= -szy + ay(y)*2*szy;
func fy= -szy + (y)*2*szy;
//func fy= -szy + 4*(y-.5)*((y-.5)*(y-.5)+0*(y-.5)*(y-.5))*2*szy*szy*szy;
x=1;
cout<<" mnx mny "<<nx<<" "<<ny<<" "<<fx<<endl;
// potential function
func real barrier( real xx, real yy)
{
//return vcenter*(vx*xx*xx/szx/szx+vy*yy*yy/szy/szy);
//return -(10.9*vcenter)+vcenter/(1.0e-1+sqrt(vx*xx*xx/szx/szx+vy*yy*yy/szy/szy));
return 1*(vcenter)+vcenter*log(1.0e-1+sqrt(vx*xx*xx/szx/szx+vy*yy*yy/szy/szy));
//return 1000;
//cout<<xx<<" "<<sz2<<endl;
if ((xx<=sz2)&&(xx>= -sz2))
if ((yy<=sz2)&&(yy>= -sz2)) return vcenter;
return 0;
}
mesh Th = square(nx,ny,[fx,fy]);
fespace Vh(Th, Pk); // local finite element space
Vh v= barrier(x,y);
real cut=.1;
real tol=.5;
Th = adaptmesh(Th, 1000.0*pow(dx(v)*dx(v)+dy(v)*dy(v),3), splitpbedge=1, abserror=0, nbvx=1e4, cutoff=cut, err=tol, inquire=0, ratio=.01, hmin=.005, hmax=1.5);
medit("grid",Th);
// make everything with the new mesh
v= barrier(x,y);
medit("potential",Th,v);
// not needednow
Vh xx=x;
Vh yy=y;
fespace Wh(Th, Pk); // local finite element space
Wh psiz=1.0/(szx*szy);
Wh psi=psiz,vpsi,wpsi=psiz;
//////////////////////////////////////////////////
macro SCHM(v,w)(((Grad(w)'*Grad(v) ))) // EOM
int lowers=0;
// there has to be a better structure for this?
real[int,int] leigens(psi[].n,nstates);
real[int] energies(nstates);
// project out lower states from current solution
func real[int]proout(real[int] inpsi)
{
real[int] xxx=inpsi;
real[int] xxxf=inpsi;
for(int j=0; j<lowers; ++j)
{
for(int l=0; l<xxxf.n; ++l) { xxx[l]= leigens(l,j); }
// WTF when does it tark areeays anw ehen fespace
Wh fu; fu[]=xxxf;
Wh sht; sht[]=xxx;
real nu=int2d(Th)(sht*fu);
// this seems to be "1" now but still not sure code is right.
real de=int2d(Th)(sht*sht); // should just do once
real f=nu/de;
cout<<" j="<<j<<" f="<<f<<" nu="<<nu<<" de="<<de<<endl; cout.flush;
for(int l=0; l<psi[].n; ++l) { xxxf[l]= xxxf[l]-f*xxx[l]; }
} // j
return xxxf;
} // inpsi
for(int i=0; i<itermax; ++i)
{
// complex time stepping or filtering out the
// higher eigenvalues, almost opposite of the "power method"
// for eigenvalue finding
varf schrod( psi,vpsi,solver=CG,init=1) =
int2d(Th) ( SCHM(psi,vpsi))
+int2d(Th)( v*psi*vpsi)
+int2d(Th)( wpsi*vpsi/dt )
+on(1,2,3,4,psi=0)
;
matrix A=schrod(Wh,Wh);
real[int] rhs=schrod(0,Wh);
//psi[]=A^-1*rhs;
//set(B,solver=GMRES);
//set(A,solver=GMRES);
//set(A,solver=LU);
set(A,solver=CG);
//set(B,solver=CG,eps=.1);
//set(B,solver=CG);
//set(B,solver=CG);
//set(B,solver=LU); // nonseuqre doh
//psi[]=B^-1*rhsa;
psi[]=A^-1*rhs;
psi[]=proout(psi[]);
if (false) {
for(int jj=0; jj<lowers; ++jj)
{
// this may be ok for dot product not integral-
Wh test;
for(int k=0; k<psi[].n; ++k)
{test[][k]= leigens(k,jj); }
real nu=int2d(Th)(psi*test);
real du=int2d(Th)(test*test);
real com=nu/du;
cout<<"wtf " <<jj<<" nu "<<nu<<" du "<<du<<" eta "<<com <<endl;
for(int k=0; k<psi[].n; ++k) psi[][k]=psi[][k]-com*leigens(k,jj);
} // jj
} // false
cout.flush;
//cout<<"done solving "<<endl; cout.flush;
real pnorm =real(sqrt(int2d(Th)(psi*conj(psi))));
psi=psi/pnorm;
real l2 =real(sqrt(int2d(Th)((wpsi-psi)*conj(wpsi-psi))));
Wh H= -dxx(psi)-dyy(psi)+v*psi;
// going to all real loses some decay lifetime etc.
real ei=int2d(Th)(H*psi); // do]t product?
// this should be one after norm...
real psi2=int2d(Th)(psi*conj(psi)); // do]t product?
Wh Perp=H-ei*psi;
real perp=sqrt(int2d(Th)(Perp*conj(Perp)));
cout<<" iter= "<<i<<" lowers="<<lowers<< " e ="<<ei<<" perp="<<perp<< " psi2="<<psi2<< " pnorm="<<pnorm<< " l2="<<l2<<endl;
wpsi=psi;
if (l2<l2dun)
{
real[int] xxx=psi[];
for(int wtf=0; wtf<psi[].n; ++wtf)
{ leigens(wtf,lowers)=xxx[wtf]; }
energies(lowers)=ei;
++lowers;
wpsi=psiz;
/////////////////////////////
// project out...
///////////////////////////////
wpsi[]=proout(wpsi[]);
psi=wpsi;
} // l2 low enougth
macro myplot()cmm = "Global solution at iteration " + i, fill = 1, value = 1//
// ffglut gets sigabort after about 7 iters...
//if ((i%3)==0){ plot(amag,wait=1,myplot); }
if ((lowers!=0)||(l2<l2view)){
//NDIS(amag, " psi ");
//NDIS(psi, " psi ");
if ((l2<l2view)||((i%10)==9)){
for(int ii=0; ii<lowers; ++ii)
cout<<" energy " <<ii<<" "<<energies(ii)<<endl;
cout.flush;
ifstream file("ioctl.txt");
int ok=0;
file>> ok;
if (ok!=0) { medit("psi",Th,psi); }
}
}
} // i
[schpro2.edp|attachment](upload://jf70tVOf1C4KgnI7Idm3ze7FTy3.edp) (9.0 KB)
This seems to be the opposite