//  schroedinger.edp
//
//  Discussion:
//
//    Omega is the disk of radius R centered at (0,0).
//
//    PDE: i ut + alpha i - uxx - uyy + lambda * nu ( u ) * u = f on Omega
//    IC:  u(x,y,0) = u0 in Omega
//    BC:  u(x,y,t) = 0 on dOmega, for t > 0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 August 2015
//
//  Author:
//
//    Georges Sadaka
//
cout << "\n";
cout << "schroedinger:\n";
cout << "  FreeFem++ version\n";
cout << "  Set up and solve the Schroedinger equation.\n";

verbosity = 0;
//
//  Define the border of the region.
//
real R = 10.0;
border C ( t = 0.0, 2.0 * pi ) { x = R * cos(t); y = R*sin(t); label = 1; };
//
//  Mesh the region.
//
real Dx = 0.2;
int n = floor ( 2.0 * pi * R / Dx );
mesh Th = buildmesh ( C ( n ) );
//
//  Plot the mesh.
//
plot ( Th, wait = false, ps = "schroedinger_mesh.ps" );
//
//  Define the finite element space.
//
fespace Vh(Th,P2);

real dt=0.01;
real Tf=5.0;
real lambda=-1.0;
real p=3.0;
real alpha=0.0;

Vh <complex> uh;
Vh <complex> vh;
Vh <complex> uh0=exp(-x^2-y^2-5.*1i*x);
Vh <complex> uhk=uh0;
Vh <complex> TNL;
Vh <complex> B;

varf a(u,v) = int2d(Th)(u*v*1i/dt + u*v*1i*alpha/2. + (dx(u)*dx(v) + dy(u)*dy(v))/2.) + on(1,u=0);

matrix<complex> A = a(Vh,Vh);

varf b(u,v) = 
    int2d(Th) ( 
      uh0 * v * 1i / dt 
    - uh0 * v * 1i * alpha / 2.0 
    - ( dx(uh0)*dx(v) + dy(uh0)*dy(v) )/2.0 
    - lambda*TNL*(uhk+uh0)*v/2.0 ) 
  + on ( 1, u = 0.0 );

Vh ABSU;
int kk=0;
real[int] NORML2(floor(Tf/dt)+1);

for ( real t = 0.0; t <= Tf; t = t + dt )
{
  TNL = abs(uh0)^(p-1);

  for ( int i = 0; i < 2; i++ )
  {
    B[] = b(0,Vh);
    set ( A, solver = sparsesolver );
    uhk[] = A^-1 * B[];
  }

  uh0[]=uhk[];
  ABSU=abs(uh0);
  NORML2[kk]=sqrt(int2d(Th)(abs(uh0)^2));        

  if ( ! ( kk % 50 ) )
  {
    plot ( ABSU, fill = true, value = true, dim = 2, 
      ps = "schroedinger"+kk+".ps", cmm="t="+t+"  ;||u||_L^2="+NORML2[kk] );
    {
      ofstream gnufile ( "u_norm.gnu");
      for ( int i = 0; i <= kk; i++ )
      {
        gnufile << i*dt << " " << NORML2(i) << endl;
      }
    }
    exec ( 
      "echo 'plot \"u_norm.gnu\" w lp \
       pause 5 \
       quit' | gnuplot" );
  }
  kk++;
}
//
//  Terminate.
//
cout << "\n";
cout << "schroedinger:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );