Yes!!
I think I did it…Now the code is running, and it’s something new!
(It will finish in several hours, I´ll tell you if it definitely works)
Here is my code:
from dolfin import *
import matplotlib.pyplot as plt
set_log_level(LogLevel.ERROR)
mesh_NS = Mesh()
with XDMFFile("mesh_NS.xdmf") as infile:
infile.read(mesh_NS)
mvc = MeshValueCollection("size_t", mesh_NS, 2)
with XDMFFile("mf_NS.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf_NS = cpp.mesh.MeshFunctionSizet(mesh_NS, mvc)
T =70 # final time
num_steps =8000 # number of time steps
dt = T / num_steps # time step size
mu = 0.000018 # dynamic viscosity
rho = 1.20 # density
# Define function spaces
V = VectorFunctionSpace(mesh_NS, 'P', 2)#para la velocidad
Q = FunctionSpace(mesh_NS, 'P', 1)#para la presión
#Boundaries
#Los tenemos definidos ya en los grupos físicos de la malla:
#11:'aerogenerador'
#12:'fronteras superior e inferior'
#13:'entrada flujo'
#14:'salida flujo'
# 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)
u_1,u_2=split(u)
u_1=Function(V)
u_2=Function(V)
u_11=Function(V)
u_12=Function(V)
u_22=Function(V)
u_11=inner(u_1,u_1)
u_12=inner(u_1,u_2)
u_22=inner(u_2,u_2)
# Define boundary conditions
U_x=0.044
bcu_inflow = DirichletBC(V, Constant((U_x, 0)),mf_NS, 13)
bcu_cylinder = DirichletBC(V, Constant((0, 0)),mf_NS, 11)
bcu_walls = DirichletBC(V.sub(1), Constant(0),mf_NS, 12) #symmetry conditions
bcp_outflow = DirichletBC(Q, Constant(0),mf_NS, 14)
bcu = [bcu_inflow,bcu_cylinder,bcu_walls]
bcp = [bcp_outflow]
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh_NS)
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('NS_and_TL/velocitySC_200.xdmf')
xdmffile_p = XDMFFile('NS_and_TL/pressureSC_200.xdmf')
xdmffile_TL = XDMFFile('NS_and_TL/Lighthill_tensor_200.xdmf')
# Create time series (for use in acoustic .py)
timeseries_u = TimeSeries('NS_and_TL/velocitySC_200_series')
timeseries_p = TimeSeries('NS_and_TL/pressureSC_200_series')
timeseries_TL = TimeSeries('NS_and_TL/Lighthill_tensor_200_series')
# Create progress bar
progress = Progress('Time-stepping', num_steps)
# 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')
#Resuelvo tensor de Lighthill
# Tensor_Lighthill
def TL(u_11,u_12,u_22):
"Return de Lighthill Tensor projected into same space as u"
u_1,u_2=split(u)
u_1=Function(V)
u_2=Function(V)
u_11=Function(V)
u_12=Function(V)
u_22=Function(V)
u_11=inner(u_1,u_1)
u_12=inner(u_1,u_2)
u_22=inner(u_2,u_2)
degree=V.ufl_element().degree()
L=VectorFunctionSpace(mesh_NS,'P',degree)
TL=project(as_vector([rho*(u_11.dx(0)+u_12.dx(1)),rho*(u_12.dx(0)+u_22.dx(1))]),L)
return TL
# Save solution to file (XDMF/HDF5) Para visualización
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
xdmffile_TL.write(TL(u_11,u_12,u_22), t)
# Save nodal values to file (para importar en otro .py)
timeseries_u.store(u_.vector(), t)
timeseries_p.store(p_.vector(), t)
timeseries_TL.store(TL(u_11,u_12,u_22).vector(), t)
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
# Update progress bar
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
print('u max:', u_.vector().get_local().max())
print('p max:', p_.vector().get_local().max())
I did some changes, but the most important, appart from the definition of TL() as a function of u_11_u_12 and u_22 is that I had to define u_11,u_12,u_22 in the beginning too, and when a I tried to write the TL to the xdmffile and timeseries, I had to write TL(u_11,u_12,u_22) (I agree with you dokken, I should avoid different names for the same function).
Many thanks to all of you!
is not equal to
=