//  driven_cavity_navier_stokes_pod.edp
//
//  Discussion:
//
//    Proper Orthogonal Decomposition for a fluid dynamics problem.
//
//    in this example the snapshots are a series of steady Navier-Stokes
//    solutions, on a 2D lid-driven cavity. The physical parameter varied
//    is the Reynolds number.
//
//    The algorithm is as follows:
//    1) compute the snapshots
//    2) compute the modes
//    3) compute the Galerkin projection on the reduced order basis
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/driven_cavity_navier_stokes_pod/driven_cavity_navier_stokes_pod.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 May 2020
//
//  Author:
//
//    Giuseppe Pitton
//
  cout << "\n";
  cout << "driven_cavity_navier_stokes_pod:\n";
  cout << "  FreeFem++ version\n";
  cout << "  Proper Orthogonal Decomposition (POD) for a fluid flow problem\n";
  cout << "  governed by the Navier-Stokes equations.\n";

  load "lapack"
//
//  Mesh the unit square.
//
  int n = 15;
  mesh Th = square ( n, n );
//
//  Use a density function to make the mesh finer near the boundaries.
//
  Th = movemesh ( Th, [ (1.0-cos(pi*x))/2.0, (1.0-cos(pi*y))/2.0 ] );
  plot ( Th, wait = true, ps = "driven_cavity_navier_stokes_pod_mesh.ps" );
//
//  Define the finite element space for velocities and pressure.
//
  fespace Vh ( Th, [ P2, P2, P1 ] );
//
//  Define three solution vectors.
//
  Vh [u,v,p];
  Vh [uu,vv,pp];
  Vh [up,vp,q];
  up[] = 0.0;
//
// Define an array of snapshot solution vectors, of size NSNSH.
//
  int nsnsh = 4;
  Vh[int] [usnsh,vsnsh,psnsh](nsnsh);
//
//  Define a vector of average solution values, 
//  and initialize all components to 0.
//
//  "uavg[]" refers to the entire array [uavg,vavg,pavg].
//
  Vh [uavg,vavg,pavg];
  uavg[] = 0.0;
//
//  Define the divergence of a velocity vector.
//
  macro div(u,v) ( dx(u) + dy(v) )
//
//  Define the gradient vector of a scalar (pressure).
//
  macro grad(u) [ dx(u), dy(u) ]
//
//  Define the weak form of the unsteady Navier Stokes equations.
//
  real Re = 1.0;
  real nu = 1.0 / Re;
  real dt = 0.1;

  problem NSunst ( [u,v,p], [uu,vv,pp], solver = sparsesolver ) =
      int2d(Th) ( nu * ( grad(u)' * grad(uu) + grad(v)' * grad(vv) ) )
    + int2d(Th) ( [up,vp]' * grad(u) * uu + [up,vp]' * grad(v) * vv )
    - int2d(Th) ( p * div(uu,vv) + pp * div(u,v))
    + int2d(Th) ( 1.0e-10 * p * pp )
    + on ( 1,2,4, u = 0.0, v = 0.0 )
    + on ( 3, u = 1.0, v = 0.0 );
//
//  Generate the snapshots, one at each Reynolds number.
//
  for ( int j = 0; j < nsnsh; j++ )
  {
//
//  Fixed point iteration for the nonlinearity.
//
    for ( int i = 0; i < 10; i++ )
    {
      NSunst;
      up[] = u[];
    }
    usnsh[j][] = u[];
    uavg[] = uavg[] + usnsh[j][];

    plot ( [usnsh[j],vsnsh[j]], cmm = "Reynolds="+Re );

    Re = Re + 50.0;
    nu = 1.0 / Re;
  }

  uavg[] = uavg[] / nsnsh;

  plot ( [uavg,vavg], value = true, wait = false, 
    cmm = "uavg", 
    ps = "driven_cavity_navier_stokes_pod_uavg.ps" );
//
//  Subtract off the average value.
//
  for ( int i = 0; i < nsnsh; i++ )
  {
    usnsh[i][] = usnsh[i][] - uavg[];
  }
//
//  Build the correlation matrix C_ij = (u_i,u_j)
//
  real[int,int] C(nsnsh,nsnsh);

  for ( int i = 0; i < nsnsh; i++ )
  {
    for ( int j = i; j < nsnsh; j++ )
    {
      C(j,i) = int2d(Th) ( usnsh[i]*usnsh[j] + vsnsh[i]*vsnsh[j] );
    }
  }
//
//  Fill in symmetric part.
//
  for ( int j = 0; j < nsnsh; j++ )
  {
    for ( int i = j; i < nsnsh; i++ )
    {
      C(j,i) = C(i,j);
    }
  }

  cout << "correlation matrix: "<< endl;
  cout << C << endl;
//
//  Create an identity matrix.
//
  real[int,int] II(nsnsh,nsnsh);

  II = 0.0; 
  for ( int i = 0; i < nsnsh; i++ )
  {
    II(i,i) = 1.0;
  }
//
//  Call the LAPACK function that computes eigenvalues and eigenvectors.
//
  real[int] ev(nsnsh);
  real[int,int] eVec(nsnsh,nsnsh);

  int k1 = dsygvd ( C, II, ev, eVec );

  cout << "Eigenvalues of the correlation matrix: " << ev << endl;
//
//  Create a matrix containing a copy of the snapshots.
//
  real[int,int] matVec(nsnsh,Vh.ndof);

  for ( int i = 0; i < nsnsh; i++ )
  {
    matVec(i,:) = usnsh[i][];  
  }
//
//  Recover NEV modes, normalize them, and store them in the usnsh[][] matrix
//
  int nev = 4;
  real normsn;

  for ( int i = 0; i < nev; i++ )
  {
    usnsh[i][] = 0.0;
    for ( int j = 0; j < nsnsh; j++ )
    {
      usnsh[i][] = usnsh[i][] + eVec(j,nsnsh-1-i) * matVec(j,:);
    }
    normsn = sqrt ( int2d(Th) ( square ( usnsh[i] ) + square ( vsnsh[i] ) ) );
    usnsh[i][] = usnsh[i][] / normsn;
  }

  for ( int i = 0; i < nev; i++ )
  {
    plot ( [usnsh[i],vsnsh[i]], wait = false, 
      cmm = "mode"+i, ps = "driven_cavity_navier_stokes_pod_mode"+i+".ps" );
  }
//
//  Build reduced-order matrices for the online phase.
//  M is the mass matrix, M_ij=(phi_i,phi_j)
//  K is the stiffness matrix, K_ij=(grad(phi_i),grad(phi_j))
//  C is the convection tensor, C_ijk=(phi_i,phi_j.grad(phi_k))
//
  cout << "Building matrices..." << endl;

  real Cijk;
  real[int,int] M(nev,nev);
  real[int,int] K(nev,nev);
  real[int,int] Clin(nev,nev^2);

  for ( int i = 0; i < nev; i++ )
  {
    cout << i + 1 << "\/" << nev << endl;
    for ( int j = 0; j < nev; j++ )
    {
      M(i,j) = int2d(Th) ( usnsh[i]*usnsh[j] + vsnsh[i]*vsnsh[j] );
      K(i,j) = int2d(Th) ( grad ( usnsh[i] )' * grad ( usnsh[j] ) 
                         + grad ( vsnsh[i] )' * grad ( vsnsh[j] ) );
      for ( int k = 0; k < nev; k++ )
      {
        Cijk = int2d(Th) ( usnsh[i] * ( [usnsh[j],vsnsh[j]]' * grad(usnsh[k]) )
                         + vsnsh[i] * ( [usnsh[j],vsnsh[j]]' * grad(vsnsh[k]) ) );
        Clin(i,j*nev+k) = Cijk;
      }
    }
  }

  real[int,int] C1(nev,nev);
  real[int,int] C2(nev,nev);
  real[int] C3(nev);
  real[int] Kav(nev);

  for ( int i = 0; i < nev; i++ )
  {
    for ( int j = 0; j < nev; j++ )
    {
      C1(i,j) = int2d(Th) ( [uavg,vavg]' * grad(usnsh[j] ) * usnsh[i]
                          + [uavg,vavg]' * grad(vsnsh[j] ) * vsnsh[i] );

      C2(i,j) = int2d(Th) ( [usnsh[j],vsnsh[j]]' * grad(uavg)*usnsh[i]
                          + [usnsh[j],vsnsh[j]]' * grad(vavg)*vsnsh[i] );
    }
    Kav(i) = int2d(Th) ( grad(usnsh[i])' * grad(uavg) 
                       + grad(vsnsh[i])' * grad(vavg) );

    C3(i) = int2d(Th) ( [uavg,vavg]' * grad(uavg) * usnsh[i]
                      + [uavg,vavg]' * grad(vavg) * vsnsh[i] );
  }
//
//  Solve the reduced system at Reynolds number RET.
//
  real Ret = 101.0;
  nu = 1.0 / Ret;
  cout << "Solving reduced system..." << endl;
  cout << "Target Reynolds number = " << Ret << endl;

  real[int] a(nev);
  real[int] ap(nev);
  real[int] bb(nev);
  real[int,int] MM(nev,nev);
  real[int,int] Minv(nev,nev);

  MM = 0.0;
  bb = 0.0;
  a = 0.0;
  ap = 0.0;

  for ( int k = 0; k < 100; k++ )
  {
    real residual = 0.0;
    for ( int i = 0; i < nev; i++ )
    {
      for ( int j = 0; j < nev; j++ )
      {
        MM(i,j) = nu * K(i,j) + C1(i,j) + C2(i,j);
        for ( int l = 0; l < nev; l++ )
        {
          MM(i,j) = MM(i,j) + Clin(i,j*nev+l) * ap(l);
        }
      }
      bb(i) = C3(i) + nu * Kav(i);
    }
//
//  Compute inverse(MM), to solve MM * a = bb.
//
    Minv = MM^-1;
    a = Minv * bb;
    for ( int ii = 0; ii < nev; ii++ )
    {
      residual = residual + ( a(ii) - ap(ii) )^2;
    }
    residual = sqrt ( residual );
    cout << residual << endl;
    ap = a;

    if ( residual < 1.0e-8 ) cout << "tolerance reached" << endl;
    if ( residual < 1.0e-8 ) break;
  }
  cout << "Computed coefficients: " << a << endl;
//
//  Reduced order solutions
//
  Vh [urebuilt,vrebuilt,prebuilt] = [0.0,0.0,0.0];

  urebuilt[] = uavg[];
  for ( int i = 0; i < nev; i++ )
  {
    urebuilt[] = urebuilt[] + a(i) * usnsh[i][];
  }

  plot ( [urebuilt,vrebuilt], wait = true,
    cmm = "Reduced order solution, Re=" + Ret );
//
//  Terminate.
//
cout << "\n";
cout << "driven_cavity_navier_stokes_pod:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );