Startpage >> Main >> POD

POD

POD for Navier-Stokes (thanks to Giuseppe Pitton, SISSA mathLab, March 2014, gpitton@sissa.it):

In this script we show how to build an offline/online computational reduction using Proper Orthogonal Decomposition (POD) in FreeFem++. As example problem, we choose the standard 2D cavity flow. It consists in solving the following steady state Navier-Stokes equations:

Find \((\mathbf{u},p)\) such that: \begin{array}{rcl} \mathbf{u}\cdot\nabla\mathbf{u}-\nu\Delta\mathbf{u}+\nabla p=\mathbf{0}\\ \nabla\cdot\mathbf{u}=0, \end{array} with \(\nu=\frac{1}{\mathrm{Re}}\), where \(\mathrm{Re}\) is the Reynolds number. We consider non-homogeneous Dirichlet boundary conditions on the upper side, homogeneous Dirichlet (no-slip) boundary conditions on the remaining sides.

For the nonlinearity we adopt a fixed-point iteration scheme, that in variational form reads:

For each step \(k\), find \(\mathbf{u}^k\in\mathbf{H}^1_0\) and \(p^k\in L^2\) such that \begin{array}{rcl} & (\mathbf{v},\mathbf{u}^{k-1}\cdot\nabla\mathbf{u}^k)+\nu\,(\nabla\mathbf{v},\nabla\mathbf{u}^k)-(\nabla\cdot\mathbf{v}, p^k)=\mathbf{0}\qquad&\forall\mathbf{v}\in\mathbf{H}^1_0\\ &(q,\nabla\cdot\mathbf{u}^k)=0 &\forall q\in L^2 \end{array} until a tolerance is reached, namely until \(\|\mathbf{u}^k-\mathbf{u}^{k-1}\|<\mathrm{tol}\).

The main goal is to rapidly solve the equations for any value of the Reynolds number inside a given parameter range \((\mathrm{Re}_a,\mathrm{Re}_b)\). For this reason, first (during the offline phase) we compute a sequence of representative solutions (called snapshots) \((\mathbf{u}_i,p_i)_{i=1}^{N_{sn}}\) for a given Reynolds number \(\mathrm{Re}_i\). Then, we build a basis set \((\mathbf{\varphi}_i,\psi_i)_{i=1}^N\) for the velocity and pressure spaces that will be used to build the Galerkin projection during the online phase:

For each step \(k\), find \((\mathbf{u}^{N,k},p^{N,k})\) such that

\begin{array}{rcl} & \sum_{j=1}^N u_j^{N,k}(\mathbf{\varphi}_i,\sum_{r=1}^N\mathbf{\varphi}_r\cdot\nabla\mathbf{\varphi}_j)+\nu\, \sum_{j=1}^N u_j^{N,k}(\nabla\mathbf{\phi}_i,\nabla\mathbf{\varphi}_j)- \sum_{j=1}^N p_j^{N,k}(\nabla\cdot\mathbf{\varphi}_i, \psi_j) =\mathbf{0} \qquad & \forall \mathbf{\varphi}_i\\ & \sum_{j=1}^N u_j^{N,k}(\psi_i,\nabla\cdot\mathbf{\varphi}_j)=0 &\forall \psi_i\in L^2 \end{array}

Since the basis \((\mathbf{\varphi}_i,\psi_i)\) has a much lower dimensionality than the (local) basis used during the offline phase, the solution of this new problem is expected to be much faster than with the standard Finite Element Method.

In this example, the reduced order basis is obtained by means of a Proper Orthogonal Decomposition technique, that consists first in forming the correlation matrix \(C\in\mathbb{R}^{N_{sn}^2}\) of the snapshots: \( C_{ij}=(\mathbf{u}_i,\mathbf{u}_j) \)

Then, the eigenpairs \((\lambda_k,\mathbf{w}_k)\) of \(C\) are computed, and used to build the basis functions: \( \mathbf{\phi}_i=\sum_{n=1}^{N_{sn}}\mathbf{w}_n\mathbf{u}_n \) The basis functions so obtained may be normalized. Other important features of the POD is that the basis functions \(\mathbf{\phi}_i\) are \(L^2\)-orthogonal and fullfil a best-approximation property in the \(L^2\) norm.

POD Solution for Re=101

download example: podNS.edp or return to main page

/*********************************************************************
 * 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
 * usage: FreeFem++ podNS.edp
 * written by Giuseppe Pitton, SISSA mathLab, March 2014
 * gpitton@sissa.it
 *********************************************************************/
load "lapack"

int n=15;    // mesh refinement
int nev=4;   // number of modes to compute
int nsnsh=4; // number of snapshots to use

real Ret=101;     // target Reynolds number (for the offline phase)

mesh Th=square(n,n);
Th=movemesh(Th,[(1-cos(pi*x))/2.,(1-cos(pi*y))/2.]);  // refine near the boundaries
fespace Vh(Th,[P2,P2,P1]);
Vh [u,v,p],[uu,vv,pp],[up,vp,q];
// snapshot space:
Vh[int] [usnsh,vsnsh,psnsh](nsnsh);
// average values:
Vh [uavg,vavg,pavg];
uavg[]=0.; up[]=0.;

macro div(u,v)(dx(u)+dy(v))//
macro grad(u)[dx(u),dy(u)]//

int i,j,k;
real Re=1,nu=1/Re,dt=.1;
real[int] viso=[0.33,0.21,0.17,0.12,0.07,0.03,-0.017,-0.0065,
-0.112,-0.16,-0.21,-0.25,-0.3,-0.35,-0.4,-0.44,-0.49,-0.54,-0.585,-0.66];

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.e-10*p*pp)
  +on(1,2,4,u=0,v=0)
  +on(3,u=1,v=0)
  ;

/*for(int i=0;i<60;i++){
  NSunst;
  up[]=u[];
}*/

for (j=0;j<nsnsh;j++){
  for (i=0;i<10;i++){ // here starts the fixed-point iteration for the nonlinearity
    NSunst;
    up[]=u[];
  }
  usnsh[j][]=u[];
  uavg[]+=usnsh[j][];
  plot([usnsh[j],vsnsh[j]],cmm="Reynolds="+Re);
  // increase the Reynolds number by 50 for each snapshot
  Re += 50;nu=1/Re;
}

uavg[]/=nsnsh;
plot([uavg,vavg],value=1,cmm="uavg",wait=0);

for (int i=0;i<nsnsh;i++)
  usnsh[i][]-=uavg[];

/* build the correlation matrix:
 * C_ij = (u_i,u_j)
 */
real[int,int] C(nsnsh,nsnsh);
for (i=0;i<nsnsh;i++)
  for (j=i;j<nsnsh;j++)
    C(j,i) = int2d(Th)(usnsh[i]*usnsh[j]+vsnsh[i]*vsnsh[j]);
for (j=0;j<nsnsh;j++)
  for (i=j;i<nsnsh;i++)
    C(j,i)=C(i,j);

cout << "correlation matrix: "<< endl;
cout << C << endl;

// identity matrx
real[int,int] II(nsnsh,nsnsh);
II=0; 
for (i=0;i<nsnsh;i++) II(i,i)=1.;

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;

real[int,int] matVec(nsnsh,Vh.ndof);
for (i=0;i<nsnsh;i++)
  matVec(i,:)=usnsh[i][];  

real normsn=1;
// recover the modes and store them in the usnsh[][] matrix
for (i=0;i<nev;i++){
  usnsh[i][]=0.;
  for (j=0;j<nsnsh;j++){
    usnsh[i][]+=eVec(j,nsnsh-1-i)*matVec(j,:); // reordering: most energetic modes first
  }
  normsn = sqrt(int2d(Th)(square(usnsh[i])+square(vsnsh[i])));
  usnsh[i][]/=normsn;
}


plot([uavg,vavg],wait=0);
for(i=0;i<nev;i++) {
  plot([usnsh[i],vsnsh[i]],cmm="mode"+i,ps="PODmode"+i+".eps",wait=1);
}

// building reduced-order matrices for the online phase
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); // tensor form of C
for (i=0;i<nev;i++){
  cout << i+1 << "\/" << nev << endl;
  for (j=0;j<nev;j++){
    /* mass M_ij=(phi_i,phi_j) */
    M(i,j)=int2d(Th)(usnsh[i]*usnsh[j]+vsnsh[i]*vsnsh[j]);
    /* stiffness K_ij=(grad(phi_i),grad(phi_j)) */
    K(i,j)=int2d(Th)(grad(usnsh[i])'*grad(usnsh[j])+grad(vsnsh[i])'*grad(vsnsh[j]));
    for (k=0;k<nev;k++){
    /* convection: C_ijk=(phi_i,phi_j.grad(phi_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),C2(nev,nev);
real[int]C3(nev),Kav(nev);
for (i=0;i<nev;i++){
  for (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 reduced system */
nu=1/Ret;
cout << "Solving reduced system..." << endl;
cout << "Target Reynolds number = " << Ret << endl;
real[int,int] MM(nev,nev),Minv(nev,nev);
real[int] bb(nev),a(nev),ap(nev);
MM=0; bb=0; a=0; ap=0;
for (k=0;k<100;k++){
real residual=0.;
for (i=0;i<nev;i++){
  for (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)+=Clin(i,j*nev+l)*ap(l);
  }
  bb(i)=C3(i)+nu*Kav(i);
}
Minv=MM^-1;
a=Minv*bb;
for(int ii=0;ii<nev;ii++)
   residual+=(a(ii)-ap(ii))^2;
residual=sqrt(residual);
cout << residual << endl;
ap=a;
if (residual<1.e-8) cout << "tolerance reached" << endl;
if (residual<1.e-8) break;
}
cout << "computed coefficients: " << a << endl;

/* reduced order solutions */
Vh [urebuilt,vrebuilt,prebuilt]=[0,0,0];

urebuilt[]=uavg[];
for (int i=0;i<nev;i++)
   urebuilt[]+=a(i)*usnsh[i][];
plot([urebuilt,vrebuilt],wait=1,cmm="Reduced order solution, Re="+Ret);

return to main page

Page last modified on April 12, 2016, at 01:41 PM