-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // driven_cavity_navier_stokes_pod.edp 2 : // 3 : // Discussion: 4 : // 5 : // Proper Orthogonal Decomposition for a fluid dynamics problem. 6 : // 7 : // in this example the snapshots are a series of steady Navier-Stokes 8 : // solutions, on a 2D lid-driven cavity. The physical parameter varied 9 : // is the Reynolds number. 10 : // 11 : // The algorithm is as follows: 12 : // 1) compute the snapshots 13 : // 2) compute the modes 14 : // 3) compute the Galerkin projection on the reduced order basis 15 : // 16 : // Location: 17 : // 18 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/driven_cavity_navier_stokes_pod/driven_cavity_navier_stokes_pod.edp 19 : // 20 : // Modified: 21 : // 22 : // 25 May 2020 23 : // 24 : // Author: 25 : // 26 : // Giuseppe Pitton 27 : // 28 : cout << "\n"; 29 : cout << "driven_cavity_navier_stokes_pod:\n"; 30 : cout << " FreeFem++ version\n"; 31 : cout << " Proper Orthogonal Decomposition (POD) for a flu ... : id flow problem\n"; 32 : cout << " governed by the Navier-Stokes equations.\n"; 33 : 34 : load "lapack" Add lapack interface ... 35 : // 36 : // Mesh the unit square. 37 : // 38 : int n = 15; 39 : mesh Th = square ( n, n ); 40 : // 41 : // Use a density function to make the mesh finer near the boundaries. 42 : // 43 : Th = movemesh ( Th, [ (1.0-cos(pi*x))/2.0, (1.0-cos(pi*y))/2.0 ] ); 44 : plot ( Th, wait = 1, ps = "driven_cavity_navier_stokes_pod_mesh.ps" ); 45 : // 46 : // Define the finite element space for velocities and pressure. 47 : // 48 : fespace Vh ( Th, [ P2, P2, P1 ] ); 49 : // 50 : // Define three solution vectors. 51 : // 52 : Vh [u,v,p]; 53 : Vh [uu,vv,pp]; 54 : Vh [up,vp,q]; 55 : up[] = 0.0; 56 : // 57 : // Define an array of snapshot solution vectors, of size NSNSH. 58 : // 59 : int nsnsh = 4; 60 : Vh[int] [usnsh,vsnsh,psnsh](nsnsh); 61 : // 62 : // Define a vector of average solution values, 63 : // and initialize all components to 0. 64 : // 65 : // "uavg[]" refers to the entire array [uavg,vavg,pavg]. 66 : // 67 : Vh [uavg,vavg,pavg]; 68 : uavg[] = 0.0; 69 : // 70 : // Define the divergence of a velocity vector. 71 : // 72 : macro div(u,v) ( dx(u) + dy(v) ) 73 # ) // 74 : // Define the gradient vector of a scalar (pressure). 75 : // 76 : macro grad(u) [ dx(u), dy(u) ] 77 # ) // 78 : // Define the weak form of the unsteady Navier Stokes equations. 79 : // 80 : real Re = 1.0; 81 : real nu = 1.0 / Re; 82 : real dt = 0.1; 83 : 84 : problem NSunst ( [u,v,p], [uu,vv,pp], solver = sparsesolver ) = 85 : int2d(Th) ( nu * ( grad(u) 77 @ [ dx(u), dy(u) ] 77 @ ' * grad(uu) 77 @ [ dx(uu), dy(uu) ] 77 @ + grad(v) 77 @ [ dx(v), dy(v) ] 77 @ ' * grad(vv) 77 @ [ dx(vv), dy(vv) ] 77 @ ) ) 86 : + int2d(Th) ( [up,vp]' * grad(u) 77 @ [ dx(u), dy(u) ] 77 @ * uu + [up,vp]' * grad(v) 77 @ [ dx(v), dy(v) ] 77 @ * vv ) 87 : - int2d(Th) ( p * div(uu,vv) 73 @ ( dx(uu) + dy(vv) ) 73 @ + pp * div(u,v) 73 @ ( dx(u) + dy(v) ) 73 @ ) 88 : + int2d(Th) ( 1.0e-10 * p * pp ) 89 : + on ( 1,2,4, u = 0.0, v = 0.0 ) 90 : + on ( 3, u = 1.0, v = 0.0 ); 91 : // 92 : // Generate the snapshots, one at each Reynolds number. 93 : // 94 : for ( int j = 0; j < nsnsh; j++ ) 95 : { 96 : // 97 : // Fixed point iteration for the nonlinearity. 98 : // 99 : for ( int i = 0; i < 10; i++ ) 100 : { 101 : NSunst; 102 : up[] = u[]; 103 : } 104 : usnsh[j][] = u[]; 105 : uavg[] = uavg[] + usnsh[j][]; 106 : 107 : plot ( [usnsh[j],vsnsh[j]], cmm = "Reynolds="+Re ); 108 : 109 : Re = Re + 50.0; 110 : nu = 1.0 / Re; 111 : } 112 : 113 : uavg[] = uavg[] / nsnsh; 114 : 115 : plot ( [uavg,vavg], value = 1, wait = 0, cmm="uavg", ps = "driven_cavity_navier_stokes_pod_uavg.ps" ); 116 : // 117 : // Subtract off the average value. 118 : // 119 : for ( int i = 0; i < nsnsh; i++ ) 120 : { 121 : usnsh[i][] = usnsh[i][] - uavg[]; 122 : } 123 : // 124 : // Build the correlation matrix C_ij = (u_i,u_j) 125 : // 126 : real[int,int] C(nsnsh,nsnsh); 127 : 128 : for ( int i = 0; i < nsnsh; i++ ) 129 : { 130 : for ( int j = i; j < nsnsh; j++ ) 131 : { 132 : C(j,i) = int2d(Th) ( usnsh[i]*usnsh[j] + vsnsh[i]*vsnsh[j] ); 133 : } 134 : } 135 : // 136 : // Fill in symmetric part. 137 : // 138 : for ( int j = 0; j < nsnsh; j++ ) 139 : { 140 : for ( int i = j; i < nsnsh; i++ ) 141 : { 142 : C(j,i) = C(i,j); 143 : } 144 : } 145 : 146 : cout << "correlation matrix: "<< endl; 147 : cout << C << endl; 148 : // 149 : // Create an identity matrix. 150 : // 151 : real[int,int] II(nsnsh,nsnsh); 152 : 153 : II = 0.0; 154 : for ( int i = 0; i < nsnsh; i++ ) 155 : { 156 : II(i,i) = 1.0; 157 : } 158 : // 159 : // Call the LAPACK function that computes eigenvalues and eigenvectors. 160 : // 161 : real[int] ev(nsnsh); 162 : real[int,int] eVec(nsnsh,nsnsh); 163 : 164 : int k1 = dsygvd ( C, II, ev, eVec ); 165 : 166 : cout << "Eigenvalues of the correlation matrix: " << ev << endl; 167 : // 168 : // Create a matrix containing a copy of the snapshots. 169 : // 170 : real[int,int] matVec(nsnsh,Vh.ndof); 171 : 172 : for ( int i = 0; i < nsnsh; i++ ) 173 : { 174 : matVec(i,:) = usnsh[i][]; 175 : } 176 : // 177 : // Recover NEV modes, normalize them, and store them in the usnsh[][] matrix 178 : // 179 : int nev = 4; 180 : real normsn; 181 : 182 : for ( int i = 0; i < nev; i++ ) 183 : { 184 : usnsh[i][] = 0.0; 185 : for ( int j = 0; j < nsnsh; j++ ) 186 : { 187 : usnsh[i][] = usnsh[i][] + eVec(j,nsnsh-1-i) * matVec(j,:); 188 : } 189 : normsn = sqrt ( int2d(Th) ( square ( usnsh[i] ) + square ( vsnsh[i] ) ) ); 190 : usnsh[i][] = usnsh[i][] / normsn; 191 : } 192 : 193 : for ( int i = 0; i < nev; i++ ) 194 : { 195 : plot ( [usnsh[i],vsnsh[i]], wait = 0, cmm="mode"+i, ps = "driven_cavity_navier_stokes_pod_mode"+i+".ps" ); 196 : } 197 : // 198 : // Build reduced-order matrices for the online phase. 199 : // M is the mass matrix, M_ij=(phi_i,phi_j) 200 : // K is the stiffness matrix, K_ij=(grad(phi_i),grad(phi_j)) 201 : // C is the convection tensor, C_ijk=(phi_i,phi_j.grad(phi_k)) 202 : // 203 : cout << "Building matrices..." << endl; 204 : 205 : real Cijk; 206 : real[int,int] M(nev,nev); 207 : real[int,int] K(nev,nev); 208 : real[int,int] Clin(nev,nev^2); 209 : 210 : for ( int i = 0; i < nev; i++ ) 211 : { 212 : cout << i + 1 << "/" << nev << endl; 213 : for ( int j = 0; j < nev; j++ ) 214 : { 215 : M(i,j) = int2d(Th) ( usnsh[i]*usnsh[j] + vsnsh[i]*vsnsh[j] ); 216 : K(i,j) = int2d(Th) ( grad ( usnsh[i] ) 77 @ [ dx( usnsh[i]), dy( usnsh[i]) ] 77 @ ' * grad ( usnsh[j] ) 77 @ [ dx( usnsh[j]), dy( usnsh[j]) ] 77 @ 217 : + grad ( vsnsh[i] ) 77 @ [ dx( vsnsh[i]), dy( vsnsh[i]) ] 77 @ ' * grad ( vsnsh[j] ) 77 @ [ dx( vsnsh[j]), dy( vsnsh[j]) ] 77 @ ); 218 : for ( int k = 0; k < nev; k++ ) 219 : { 220 : Cijk = int2d(Th) ( usnsh[i] * ( [usnsh[j],vsnsh[j]]' * grad(usnsh[k]) 77 @ [ dx(usnsh[k]), dy(usnsh[k]) ] 77 @ ) 221 : + vsnsh[i] * ( [usnsh[j],vsnsh[j]]' * grad(vsnsh[k]) 77 @ [ dx(vsnsh[k]), dy(vsnsh[k]) ] 77 @ ) ); 222 : Clin(i,j*nev+k) = Cijk; 223 : } 224 : } 225 : } 226 : 227 : real[int,int] C1(nev,nev); 228 : real[int,int] C2(nev,nev); 229 : real[int] C3(nev); 230 : real[int] Kav(nev); 231 : 232 : for ( int i = 0; i < nev; i++ ) 233 : { 234 : for ( int j = 0; j < nev; j++ ) 235 : { 236 : C1(i,j) = int2d(Th) ( [uavg,vavg]' * grad(usnsh[j] ) 77 @ [ dx(usnsh[j]), dy(usnsh[j]) ] 77 @ * usnsh[i] 237 : + [uavg,vavg]' * grad(vsnsh[j] ) 77 @ [ dx(vsnsh[j]), dy(vsnsh[j]) ] 77 @ * vsnsh[i] ); 238 : 239 : C2(i,j) = int2d(Th) ( [usnsh[j],vsnsh[j]]' * grad(uavg) 77 @ [ dx(uavg), dy(uavg) ] 77 @ *usnsh[i] 240 : + [usnsh[j],vsnsh[j]]' * grad(vavg) 77 @ [ dx(vavg), dy(vavg) ] 77 @ *vsnsh[i] ); 241 : } 242 : Kav(i) = int2d(Th) ( grad(usnsh[i]) 77 @ [ dx(usnsh[i]), dy(usnsh[i]) ] 77 @ ' * grad(uavg) 77 @ [ dx(uavg), dy(uavg) ] 77 @ 243 : + grad(vsnsh[i]) 77 @ [ dx(vsnsh[i]), dy(vsnsh[i]) ] 77 @ ' * grad(vavg) 77 @ [ dx(vavg), dy(vavg) ] 77 @ ); 244 : 245 : C3(i) = int2d(Th) ( [uavg,vavg]' * grad(uavg) 77 @ [ dx(uavg), dy(uavg) ] 77 @ * usnsh[i] 246 : + [uavg,vavg]' * grad(vavg) 77 @ [ dx(vavg), dy(vavg) ] 77 @ * vsnsh[i] ); 247 : } 248 : // 249 : // Solve the reduced system at Reynolds number RET. 250 : // 251 : real Ret = 101.0; 252 : nu = 1.0 / Ret; 253 : cout << "Solving reduced system..." << endl; 254 : cout << "Target Reynolds number = " << Ret << endl; 255 : 256 : real[int] a(nev); 257 : real[int] ap(nev); 258 : real[int] bb(nev); 259 : real[int,int] MM(nev,nev); 260 : real[int,int] Minv(nev,nev); 261 : 262 : MM = 0.0; 263 : bb = 0.0; 264 : a = 0.0; 265 : ap = 0.0; 266 : 267 : for ( int k = 0; k < 100; k++ ) 268 : { 269 : real residual = 0.0; 270 : for ( int i = 0; i < nev; i++ ) 271 : { 272 : for ( int j = 0; j < nev; j++ ) 273 : { 274 : MM(i,j) = nu * K(i,j) + C1(i,j) + C2(i,j); 275 : for ( int l = 0; l < nev; l++ ) 276 : { 277 : MM(i,j) = MM(i,j) + Clin(i,j*nev+l) * ap(l); 278 : } 279 : } 280 : bb(i) = C3(i) + nu * Kav(i); 281 : } 282 : // 283 : // Compute inverse(MM), to solve MM * a = bb. 284 : // 285 : Minv = MM^-1; 286 : a = Minv * bb; 287 : for ( int ii = 0; ii < nev; ii++ ) 288 : { 289 : residual = residual + ( a(ii) - ap(ii) )^2; 290 : } 291 : residual = sqrt ( residual ); 292 : cout << residual << endl; 293 : ap = a; 294 : 295 : if ( residual < 1.0e-8 ) cout << "tolerance reached" << endl; 296 : if ( residual < 1.0e-8 ) break; 297 : } 298 : cout << "Computed coefficients: " << a << endl; 299 : // 300 : // Reduced order solutions 301 : // 302 : Vh [urebuilt,vrebuilt,prebuilt] = [0.0,0.0,0.0]; 303 : 304 : urebuilt[] = uavg[]; 305 : for ( int i = 0; i < nev; i++ ) 306 : { 307 : urebuilt[] = urebuilt[] + a(i) * usnsh[i][]; 308 : } 309 : 310 : plot ( [urebuilt,vrebuilt],wait=1,cmm="Reduced order solution, Re="+Ret); 311 : // 312 : // Terminate. 313 : // 314 : cout << "\n"; 315 : cout << "driven_cavity_navier_stokes_pod:\n"; 316 : cout << " Normal end of execution.\n"; 317 : 318 : sizestack + 1024 =3496 ( 2472 ) driven_cavity_navier_stokes_pod: FreeFem++ version Proper Orthogonal Decomposition (POD) for a fluid flow problem governed by the Navier-Stokes equations. -- Square mesh : nb vertices =256 , nb triangles = 450 , nb boundary edges 60 -- Solve : min -540.541 max 561.619 -- Solve : min -539.974 max 562.226 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -539.975 max 562.227 -- Solve : min -10.0778 max 11.6581 -- Solve : min -10.104 max 11.694 -- Solve : min -10.1031 max 11.6928 -- Solve : min -10.1027 max 11.6925 -- Solve : min -10.1027 max 11.6925 -- Solve : min -10.1027 max 11.6925 -- Solve : min -10.1027 max 11.6925 -- Solve : min -10.1027 max 11.6925 -- Solve : min -10.1027 max 11.6925 -- Solve : min -10.1027 max 11.6925 -- Solve : min -4.88881 max 6.26988 -- Solve : min -4.90951 max 6.2969 -- Solve : min -4.9081 max 6.29418 -- Solve : min -4.90771 max 6.29404 -- Solve : min -4.90774 max 6.29416 -- Solve : min -4.90774 max 6.29417 -- Solve : min -4.90774 max 6.29416 -- Solve : min -4.90774 max 6.29416 -- Solve : min -4.90774 max 6.29416 -- Solve : min -4.90774 max 6.29416 -- Solve : min -3.15771 max 4.47102 -- Solve : min -3.17297 max 4.48817 -- Solve : min -3.172 max 4.48521 -- Solve : min -3.17209 max 4.48542 -- Solve : min -3.17219 max 4.48552 -- Solve : min -3.17218 max 4.48551 -- Solve : min -3.17218 max 4.48551 -- Solve : min -3.17218 max 4.48551 -- Solve : min -3.17218 max 4.48551 -- Solve : min -3.17218 max 4.48551 correlation matrix: 4 4 0.002432856998 0.0005809706562 -0.0009926253669 -0.002021202288 0.0005809706562 0.0002683467082 -0.0001913821968 -0.0006579351675 -0.0009926253669 -0.0001913821968 0.0004254545636 0.0007585530002 -0.002021202288 -0.0006579351675 0.0007585530002 0.001920584455 Eigenvalues of the correlation matrix: 4 3.25825724e-18 4.135022476e-06 0.0002799056849 0.004763202018 Building matrices... 1/4 2/4 3/4 4/4 Solving reduced system... Target Reynolds number = 101 0.0285929 0.000161302 1.5133e-06 2.08305e-08 2.56142e-10 tolerance reached Computed coefficients: 4 -0.01959899338 0.01248766463 0.01054572576 0.01301409543 driven_cavity_navier_stokes_pod: Normal end of execution. times: compile 0.007923s, execution 34.0702s, mpirank:0 CodeAlloc : nb ptr 4781, size :521288 mpirank: 0 Ok: Normal End