60
edits
No edit summary |
No edit summary |
||
Line 601: | Line 601: | ||
</div> | </div> | ||
<div style="clear:both"></div> | <div style="clear:both"></div> | ||
=== Incompressible Navier-Stokes equations === | |||
<div style="width: 100%;"> | |||
<div style="float:left; width: 48%"> | |||
{{Code|Define HyperElasticity problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | |||
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 = VectorElement("CG", triangle, 2) | |||
# Q = FiniteElement("CG", triangle, 1) | |||
V = FunctionSpace ('VelocityUpdate', mesh); | |||
Q = FunctionSpace ('PressureUpdate', mesh); | |||
# Define trial and test functions. From ufl file | |||
# u = TrialFunction(V) | |||
# p = TrialFunction(Q) | |||
# v = TestFunction(V) | |||
# q = TestFunction(Q) | |||
# Set parameter values. From ufl file | |||
# nu = 0.01 | |||
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. From ufl file | |||
# eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \ | |||
# nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx | |||
a1 = BilinearForm ('TentativeVelocity', V, k); | |||
# Pressure update. From ufl file | |||
# a = inner(grad(p), grad(q))*dx | |||
# L = -(1/k)*div(u1)*q*dx | |||
a2 = BilinearForm ('PressureUpdate', Q); | |||
# Velocity update | |||
# a = inner(u, v)*dx | |||
# L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx | |||
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 | |||
"Computing tentative velocity" | |||
L1 = LinearForm ('TentativeVelocity', V, k, u0, f); | |||
b1 = assemble (L1, noslip); | |||
utmp = A1 \ b1; | |||
u1 = Function ('u1', V, utmp); | |||
# Pressure correction | |||
"Computing 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 | |||
"Computing 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)); | |||
save (u1, sprintf ("u_%3.3d", i)); | |||
# Move to next time step | |||
u0 = u1; | |||
t += dt | |||
end | |||
</syntaxhighlight>}} | |||
</div> | |||
<div style="float:right; width: 48%"> | |||
{{Code|Define NS problem with fenics python|<syntaxhighlight lang="python" style="font-size:13px"> | |||
from dolfin import * | |||
# Load mesh from file | |||
mesh = Mesh("lshape.xml.gz") | |||
# Define function spaces (P2-P1) | |||
V = VectorFunctionSpace(mesh, "CG", 2) | |||
Q = FunctionSpace(mesh, "CG", 1) | |||
# Define trial and test functions | |||
u = TrialFunction(V) | |||
p = TrialFunction(Q) | |||
v = TestFunction(V) | |||
q = TestFunction(Q) | |||
# Set parameter values | |||
dt = 0.01 | |||
T = 3 | |||
nu = 0.01 | |||
# Define time-dependent pressure boundary condition | |||
p_in = Expression("sin(3.0*t)", t=0.0) | |||
# Define boundary conditions | |||
noslip = DirichletBC(V, (0, 0), | |||
"on_boundary && \ | |||
(x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \ | |||
(x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))") | |||
inflow = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS") | |||
outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS") | |||
bcu = [noslip] | |||
bcp = [inflow, outflow] | |||
# Create functions | |||
u0 = Function(V) | |||
u1 = Function(V) | |||
p1 = Function(Q) | |||
# Define coefficients | |||
k = Constant(dt) | |||
f = Constant((0, 0)) | |||
# Tentative velocity step | |||
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \ | |||
nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx | |||
a1 = lhs(F1) | |||
L1 = rhs(F1) | |||
# Pressure update | |||
a2 = inner(grad(p), grad(q))*dx | |||
L2 = -(1/k)*div(u1)*q*dx | |||
# Velocity update | |||
a3 = inner(u, v)*dx | |||
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx | |||
# Assemble matrices | |||
A1 = assemble(a1) | |||
A2 = assemble(a2) | |||
A3 = assemble(a3) | |||
# Use amg preconditioner if available | |||
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default" | |||
# Create files for storing solution | |||
ufile = File("results/velocity.pvd") | |||
pfile = File("results/pressure.pvd") | |||
# Time-stepping | |||
t = dt | |||
while t < T + DOLFIN_EPS: | |||
# Update pressure boundary condition | |||
p_in.t = t | |||
# Compute tentative velocity step | |||
begin("Computing tentative velocity") | |||
b1 = assemble(L1) | |||
[bc.apply(A1, b1) for bc in bcu] | |||
solve(A1, u1.vector(), b1, "gmres", "default") | |||
end() | |||
# Pressure correction | |||
begin("Computing pressure correction") | |||
b2 = assemble(L2) | |||
[bc.apply(A2, b2) for bc in bcp] | |||
solve(A2, p1.vector(), b2, "gmres", prec) | |||
end() | |||
# Velocity correction | |||
begin("Computing velocity correction") | |||
b3 = assemble(L3) | |||
[bc.apply(A3, b3) for bc in bcu] | |||
solve(A3, u1.vector(), b3, "gmres", "default") | |||
end() | |||
# Plot solution | |||
plot(p1, title="Pressure", rescale=True) | |||
plot(u1, title="Velocity", rescale=True) | |||
# Save to file | |||
ufile << u1 | |||
pfile << p1 | |||
# Move to next time step | |||
u0.assign(u1) | |||
t += dt | |||
print "t =", t | |||
# Hold plot | |||
interactive() | |||
</syntaxhighlight>}} | |||
</div> | |||
</div> | |||
<div style="clear:both"></div> | |||
[[Category:OctaveForge]] | [[Category:OctaveForge]] | ||
[[Category:Packages]] | [[Category:Packages]] |
edits