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

Download this file

UNS_with_L2_penalization.m    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 <gedeone-octave@users.sourceforge.net>
##
## 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 <http://www.gnu.org/licenses/>.
# 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)