## [945c19]: devel / example / Ficticious_Domain / Unsteady / L2_penalization / UNS_with_L2_penalization.m  Maximize  Restore  History

### 129 lines (105 with data), 3.6 kB

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128``` ```## Copyright (C) 2013 Marco Vassallo ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software ## Foundation; either version 3 of the License, or (at your option) any later ## version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public License along with ## this program; if not, see . # In this example we solve the 2D unsteady NS equation for a flow around a square # cylinder in a channel. In this example, the no-slip BC are computed using # a L2 penalization technique. # We use the differential Chorin-Temam algorithm for each time step. # The solution shows the typical von-karman street pattern after a transition # period. pkg load fem-fenics msh x = linspace (0, 10, 80*4+1); y = linspace (-2, 2, 80+1); msho = msh2m_structured_mesh (x, y, 1, 1:4); mesh = Mesh (msho); plot (mesh); import_ufl_Problem ("TentativeVelocity"); import_ufl_Problem ("VelocityUpdate"); import_ufl_Problem ("PressureUpdate"); TH1 = FunctionSpace ('TentativeVelocity', mesh); TH2 = FunctionSpace ('PressureUpdate', mesh); bc1 = DirichletBC (TH1, @(x, y, z, n) [0, 0], [1, 3]); bc2 = DirichletBC (TH1, @(x, y) [1.5*4*(2 - y)*(2 + y)/16 , 0], 4); bc3 = DirichletBC (TH2, @(x, y) 0, 2); bcu = {bc1, bc2}; bcp = bc3; ## Set parameter values dt = 0.1; T = 50.; u0 = Expression ('u0', @(x, y) [0; 0]); p0 = Expression ('p0', @(x, y) 0); k = Constant ('k', dt); f = Constant ('f', [0; 0]); nu_1 = 5.e-3; nu_0 = 5.e-3; r = 0.25; #ficticious domain ## cylinder # dom = @(x,y) (((x - 3).^2 + (y).^2) < (r).^2); # square dom = @(x, y) (x <= 3+r)*(x >= 3-r)*(y >= -r)*(y <= r); nu = Expression ('nu', @(x, y) dom(x, y) * nu_0 + ... (1 - dom(x, y)) * nu_1); sig_1 = 0; sig_0 = 1e8; sig = Expression ('sig',@(x, y) dom(x, y) * sig_0 + ... (1 - dom(x, y)) * sig_1); a2 = BilinearForm ('PressureUpdate', TH2, TH2); A2 = assemble (a2, bcp); a3 = BilinearForm ('VelocityUpdate', TH1, TH1); A3 = assemble(a3); ## Time-stepping t = dt; i = 0; pold = zeros (26001, 1); while t < T # Compute tentative velocity step a1 = BilinearForm ('TentativeVelocity', TH1, TH1, k, u0, sig, nu); L1 = LinearForm ('TentativeVelocity', TH1, f, p0, k, u0); [A1, b1] = assemble_system (a1, L1, bcu{:}); utmp = A1 \ b1; u1 = Function ('u1', TH1, utmp); # Pressure correction L2 = LinearForm ('PressureUpdate', TH2, k, u1); b2 = assemble (L2, bcp); ptmp = A2 \ b2; p1 = Function ('p1', TH2, ptmp); # Velocity correction L3 = LinearForm ('VelocityUpdate', TH1, k, u1, p1); b3 = assemble (L3); ut = A3 \ b3; u1 = Function ('u0', TH1, ut); #Update p p = ptmp + pold; p0 = Function ('p0', TH2, p); # Save to file save (p0, sprintf ("p_%3.3d", ++i)); delete (sprintf ("p_%3.3d.pvd", i)); save (u1, sprintf ("u_%3.3d", i)); delete (sprintf ("u_%3.3d.pvd", i)); # Move to next time step pold = p; u0 = u1; t += dt; end plot (u0); norm_err = 0; for i = 1:size(msho.p, 2) x = msho.p (1, i); y = msho.p (2, i); if (dom (x, y) == true) err_L2 = feval (u0, [x, y]); norm_err += err_L2(1).^2 + err_L2(2).^2; endif endfor error_on_bc = sqrt (norm_err) ```