-- FreeFem++ v4.6 (Thu Apr  2 15:47:38 CEST 2020 - git v4.6)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  schroedinger.edp
    2 : //
    3 : //  Discussion:
    4 : //
    5 : //    Omega is the disk of radius R centered at (0,0).
    6 : //
    7 : //    PDE: i ut + alpha i - uxx - uyy + lambda * nu ( u ) * u = f on Omega
    8 : //    IC:  u(x,y,0) = u0 in Omega
    9 : //    BC:  u(x,y,t) = 0 on dOmega, for t > 0
   10 : //
   11 : //  Licensing:
   12 : //
   13 : //    This code is distributed under the GNU LGPL license.
   14 : //
   15 : //  Modified:
   16 : //
   17 : //    20 August 2015
   18 : //
   19 : //  Author:
   20 : //
   21 : //    Georges Sadaka
   22 : //
   23 : cout << "\n";
   24 : cout << "schroedinger:\n";
   25 : cout << "  FreeFem++ version\n";
   26 : cout << "  Set up and solve the Schroedinger equation.\n";
   27 : 
   28 : verbosity = 0;
   29 : //
   30 : //  Define the border of the region.
   31 : //
   32 : real R = 10.0;
   33 : border C ( t = 0.0, 2.0 * pi ) { x = R * cos(t); y = R*sin(t); label = 1; };
   34 : //
   35 : //  Mesh the region.
   36 : //
   37 : real Dx = 0.2;
   38 : int n = floor ( 2.0 * pi * R / Dx );
   39 : mesh Th = buildmesh ( C ( n ) );
   40 : //
   41 : //  Plot the mesh.
   42 : //
   43 : plot ( Th, wait = 0, ps = "schroedinger_mesh.ps" );
   44 : //
   45 : //  Define the finite element space.
   46 : //
   47 : fespace Vh(Th,P2);
   48 : 
   49 : real dt=0.01;
   50 : real Tf=5.0;
   51 : real lambda=-1.0;
   52 : real p=3.0;
   53 : real alpha=0.0;
   54 : 
   55 : Vh <complex> uh;
   56 : Vh <complex> vh;
   57 : Vh <complex> uh0=exp(-x^2-y^2-5.*1i*x);
   58 : Vh <complex> uhk=uh0;
   59 : Vh <complex> TNL;
   60 : Vh <complex> B;
   61 : 
   62 : 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);
   63 : 
   64 : matrix<complex> A = a(Vh,Vh);
   65 : 
   66 : varf b(u,v) = 
   67 :     int2d(Th) ( 
   68 :       uh0 * v * 1i / dt 
   69 :     - uh0 * v * 1i * alpha / 2.0 
   70 :     - ( dx(uh0)*dx(v) + dy(uh0)*dy(v) )/2.0 
   71 :     - lambda*TNL*(uhk+uh0)*v/2.0 ) 
   72 :   + on ( 1, u = 0.0 );
   73 : 
   74 : Vh ABSU;
   75 : int kk=0;
   76 : real[int] NORML2(floor(Tf/dt)+1);
   77 : 
   78 : for ( real t = 0.0; t <= Tf; t = t + dt )
   79 : {
   80 :   TNL = abs(uh0)^(p-1);
   81 : 
   82 :   for ( int i = 0; i < 2; i++ )
   83 :   {
   84 :     B[] = b(0,Vh);
   85 :     set ( A, solver = sparsesolver );
   86 :     uhk[] = A^-1 * B[];
   87 :   }
   88 : 
   89 :   uh0[]=uhk[];
   90 :   ABSU=abs(uh0);
   91 :   NORML2[kk]=sqrt(int2d(Th)(abs(uh0)^2));        
   92 : 
   93 :   if ( ! ( kk % 50 ) )
   94 :   {
   95 :     plot ( ABSU, cmm="t="+t+"  ;||u||_L^2="+NORML2[kk], fill=true, value=true, dim=2, ps="schroedinger"+kk+".ps" );
   96 :     {
   97 :       ofstream gnufile ( "u_norm.gnu");
   98 :       for (int i=0;i<=kk;i++)
   99 :       {
  100 :         gnufile << i*dt << " " << NORML2(i) << endl;
  101 :       }
  102 :     }
  103 :     exec ( 
  104 :       "echo 'plot \"u_norm.gnu\" w lp \n       pause 5 \n       quit' | gnuplot" );
  107 :   }
  108 :   kk++;
  109 : }
  110 : //
  111 : //  Terminate.
  112 : //
  113 : cout << "\n";
  114 : cout << "schroedinger:\n";
  115 : cout << "  Normal end of execution.\n";
  116 : 
  117 :  sizestack + 1024 =5848  ( 4824 )


schroedinger:
  FreeFem++ version
  Set up and solve the Schroedinger equation.

schroedinger:
  Normal end of execution.