[d03d99]: inst / example / NavierStokes / NS.m Maximize Restore History

Download this file

NS.m    103 lines (78 with data), 2.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
pkg load fem-fenics msh
import_ufl_Problem ("TentativeVelocity");
import_ufl_Problem ("VelocityUpdate");
import_ufl_Problem ("PressureUpdate");
# We can either load the mesh from the file as in Dolfin but we can also use the msh pkg to generate the L-shape domain
name = [tmpnam ".geo"];
fid = fopen (name, "w");
fputs (fid,"Point (1) = {0, 0, 0, 0.1};\n");
fputs (fid,"Point (2) = {1, 0, 0, 0.1};\n");
fputs (fid,"Point (3) = {1, 0.5, 0, 0.1};\n");
fputs (fid,"Point (4) = {0.5, 0.5, 0, 0.1};\n");
fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
fputs (fid,"Line (1) = {5, 6};\n");
fputs (fid,"Line (2) = {2, 3};\n");
fputs (fid,"Line(3) = {6,1,2};\n");
fputs (fid,"Line(4) = {5,4,3};\n");
fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
fputs (fid,"Plane Surface(8) = {7};\n");
fclose (fid);
msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
unlink (canonicalize_file_name (name));
mesh = Mesh (msho);
# Define function spaces (P2-P1). From ufl file
V = FunctionSpace ('VelocityUpdate', mesh);
Q = FunctionSpace ('PressureUpdate', mesh);
# Set parameter values
dt = 0.01;
T = 3.;
# Define boundary conditions
noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
outflow = DirichletBC (Q, @(x,y) 0, 2);
# Create functions
u0 = Expression ('u0', @(x,y) [0; 0]);
# Define coefficients
k = Constant ('k', dt);
f = Constant ('f', [0; 0]);
# Tentative velocity step.
a1 = BilinearForm ('TentativeVelocity', V, k);
# Pressure update.
a2 = BilinearForm ('PressureUpdate', Q);
# Velocity update
a3 = BilinearForm ('VelocityUpdate', V);
# Assemble matrices
A1 = assemble (a1, noslip);
A3 = assemble (a3, noslip);
# Time-stepping
t = dt; i = 0;
while t < T
# Update pressure boundary condition
inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
# Compute tentative velocity step
L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
b1 = assemble (L1, noslip);
utmp = A1 \ b1;
u1 = Function ('u1', V, utmp);
# Pressure correction
L2 = LinearForm ('PressureUpdate', Q, u1, k);
[A2, b2] = assemble_system (a2, L2, inflow, outflow);
ptmp = A2 \ b2;
p1 = Function ('p1', Q, ptmp);
# Velocity correction
L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
b3 = assemble (L3, noslip);
ut = A3 \ b3;
u1 = Function ('u0', V, ut);
# Plot solution
plot (p1);
plot (u1);
# Save to file
save (p1, 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
u0 = u1;
t += dt
end