-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae)
   file : schroedinger.edp
 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 MIT 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 = false, 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, fill = true, value = true, dim = 2, wait = true,
   96 :       ps = "schroedinger"+kk+".ps", cmm="t="+t+"  ;||u||_L^2="+NORML2[kk] );
   97 :   }
   98 :   kk++;
   99 : 
  100 : //  if ( ! ( kk % 50 ) )
  101 : //  {
  102 : //    plot ( ABSU, fill = true, value = true, dim = 2, 
  103 : //      ps = "schroedinger"+kk+".ps", cmm="t="+t+"  ;||u||_L^2="+NORML2[kk] );
  104 : //    {
  105 : //      ofstream gnufile ( "u_norm.gnu");
  106 : //      for ( int i = 0; i <= kk; i++ )
  107 : //      {
  108 : //        gnufile << i*dt << " " << NORML2(i) << endl;
  109 : //      }
  110 : //    }
  111 : //    exec ( 
  112 : //      "echo 'plot \"u_norm.gnu\" w lp \
  113 : //       pause 5 \
  114 : //       quit' | gnuplot" );
  115 : //}
  116 : //kk++;
  117 : }
  118 : //
  119 : //  Terminate.
  120 : //
  121 : cout << "\n";
  122 : cout << "schroedinger():\n";
  123 : cout << "  Normal end of execution.\n";
  124 : 
  125 : exit ( 0 );
  126 :  sizestack + 1024 =5568  ( 4544 )