ft08_navier_stokes_cylinder.py (Cylinder is shaped like a bullet)
tesisatchi opened this issue · 0 comments
Greetings, its my first time in GitHub. I am told that I could get some help over projects that i cant figure out my mistakes. The thing is a project has been asigned to me which is creating a code that solves flow around a cylinder problem. In this case the cylinder is shaped just like a bullet, a half circle plus a half rectangle combined. I have copied the code from fenicsproject.org and modified it for my use, which is shared below.
"""
FEniCS tutorial demo program: Incompressible Navier-Stokes equations
for flow around a cylinder using the Incremental Pressure Correction
Scheme (IPCS).
u' + u . nabla(u)) - div(sigma(u, p)) = f
div(u) = 0
"""
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
T = 5.0 # final time
num_steps = 5000 # number of time steps
dt = T / num_steps # time step size
mu = 0.001 # dynamic viscosity
rho = 1 # density
# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05) + Rectangle(Point(0.2,0.15), Point(0.3,0.25))
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'
# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')
# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Define functions for solutions at previous and current time steps u_n = Function(V) u_ = Function(V) p_n = Function(Q) p_ = Function(Q)`
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx
# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')
# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')
# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh
# Create progress bar
# progress = Progress('Time-stepping')
progress = Progress('Time-stepping', num_steps)
# set_log_level(PROGRESS)
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')
# Plot solution
# plot(u_, title='Velocity')
# plot(p_, title='Pressure')
# Save solution to file (XDMF/HDF5)
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
# Save nodal values to file timeseries_u.store(u_.vector(), t) timeseries_p.store(p_.vector(), t)`
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
` # Update progress bar`
` # progress.update(t / T)`
` set_log_level(LogLevel.PROGRESS)`
` progress += 1`
` set_log_level(LogLevel.ERROR)`
` print('u max:', u_.vector().get_local().max())`
# Hold plot
# interactive()
When I run this code it converges to this:
u max: 6.59379743087e+19
u max: 3.72600712496e+38
u max: 1.20155908645e+76
and then this error pops up:
RuntimeError Traceback (most recent call last)
in
131 b2 = assemble(L2)
132 [bc.apply(b2) for bc in bcp]
--> 133 solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
134
135 # Step 3: Velocity correction step
/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
225 raise RuntimeError("Not expecting keyword arguments when solving linear algebra problem.")
226
--> 227 return dolfin.la.solver.solve(*args)
228
229
/usr/local/lib/python3.6/dist-packages/dolfin/la/solver.py in solve(A, x, b, method, preconditioner)
70 """
71
---> 72 return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
I do not know what this error means, I do not know how to solve this equation and I do not know what have I done wrong. Side note, I am just a beginner in this coding environment so any positive/negative feedback is welcome.
If anyone who understands the problem please do not hesitate to help.
Sorry if i couldnt adress "insert code" operator correctly.