Hi, I’m the new to use FEniCS and working on visualizing the viscoelastic fluid flow through a cylinder.
Whenever I try to run this code, I got -nan values as the Newton residual values… Can you please check and help me to find the error of this code?
Sometimes it works for very very slow inflow velocity condition (0.00000001 * ~~) during random time step (not for whole time step), after that some time step newton residual value does not converge anymore… Furthermore, I checked the data from this specific condition with Paraview but it does not have any meaningful data…
import fenics as fe
import numpy as np
from mpi4py import MPI
import sys, os, time
rank = MPI.COMM_WORLD.Get_rank()
tag = 'A01_f1_test'
tag_dir = os.path.join('data', tag)
if rank == 0:
fe.set_log_level(20)
print("Tag: ", tag)
os.makedirs(tag_dir, exist_ok = True)
else:
fe.set_log_level(30)
sys.stdout.flush()
start = time.time()
mu = 1
rho = 1
num_cycles = 5
steps_per_cycle = 50
f = 1
dt = 1 / (f * steps_per_cycle)
num_steps = steps_per_cycle * num_cycles
n_FDI = 2
G = fe.Constant(1)
eta = fe.Constant(1)
rel = eta / G
# Define a function to save u, p , c and pressure gradient
def save_sol(c_, u_, p_, t):
strs = fe.project(stress(u_, c_), fe.TensorFunctionSpace(mesh, 'CG', 2), solver_type = 'cg', preconditioner_type = 'hypre_amg')
dp = fe.project(fe.nabla_grad(p_), fe.VectorFunctionSpace(mesh, 'CG', 1), solver_type = 'cg', preconditioner_type = 'hypre_amg')
if t == 0:
flag = False
else:
flag = True
# Save solutions to XDMF files
xdmf_u.write_checkpoint(u_, "u", t, fe.XDMFFile.Encoding.HDF5, flag)
xdmf_p.write_checkpoint(p_, "p", t, fe.XDMFFile.Encoding.HDF5, flag)
xdmf_s.write_checkpoint(strs, "s", t, fe.XDMFFile.Encoding.HDF5, flag)
xdmf_dp.write_checkpoint(dp, "dp", t, fe.XDMFFile.Encoding.HDF5, flag)
# Backward Euler
def der_be(U, Uprev, dT):
return (U - Uprev) / dT
def der_be_UCM(C, Cprev, dT):
return (C - Cprev) / dT
# Create mesh
mesh = fe.RectangleMesh(fe.Point(-0.25,-0.025), fe.Point(0.25,0.025), 10, 10, "right/left")
# Define boundaries
inflow = "near(x[0], -0.25)"
outflow = "near(x[0], 0.25)"
walls = "near(x[1], -0.025) || near(x[1], 0.025)"
# Define Functionspaces with mixed elements
P1 = fe.FiniteElement('CG', fe.triangle, degree = 2)
P2 = fe.VectorElement('CG', fe.triangle, degree = 2)
P3 = fe.TensorElement('DG', fe.triangle, degree = 2)
elem = fe.MixedElement([P1, P2, P3])
K = fe.FunctionSpace(mesh, elem)
U = fe.VectorFunctionSpace(mesh, 'CG', 2)
N = fe.TensorFunctionSpace(mesh, 'DG', 2)
z = fe.Function(K)
z_ = fe.TrialFunction(K) # used in Jacobian formulation
p, u, c = fe.split(z)
q, v, d = fe.TestFunctions(K)
u_prev = fe.Function(U)
c_prev = fe.Function(N)
# Define expressions used in variational form
n = fe.FacetNormal(mesh)
k = fe.Constant(dt)
mu = fe.Constant(mu)
rho = fe.Constant(rho)
rel = fe.Constant(rel)
# Create XDMF files for visualization output
xdmf_u = fe.XDMFFile(os.path.join(tag_dir, tag + '_u.xdmf'))
xdmf_p = fe.XDMFFile(os.path.join(tag_dir, tag + '_p.xdmf'))
xdmf_s = fe.XDMFFile(os.path.join(tag_dir, tag + '_s.xdmf'))
xdmf_dp = fe.XDMFFile(os.path.join(tag_dir, tag + '_dp.xdmf'))
# Define symmetric gradient, stress(u)
def epsilon(u):
return fe.sym(fe.nabla_grad(u))
def stress(u, c):
return G * (c - fe.Identity(len(u)))
# Defnie boundary condition for outflow pressure
def lower_boundary_pressure_fixed_point(x, boundary):
tol = 1E-15
return (abs(x[1] + 0.025) < tol) and (abs(x[0] - 0.25) < tol)
# Start time stepping: Solve
t = 0
# Define boundary condition
inflow_profile = ('0.01', '0')
bcu_inflow = fe.DirichletBC(K.sub(1), fe.Expression(inflow_profile, degree =2), inflow)
bcu_wall = fe.DirichletBC(K.sub(1), fe.Constant((0,0)), walls)
bcu_stress = fe.DirichletBC(K.sub(2), fe.Constant([[0,0], [0,0]]), inflow)
bcu_wall_stress = fe.DirichletBC(K.sub(2), fe.Constant([[0,0], [0,0]]), walls)
outflow_profile = ('0.01 * 1600 * (0.025 * 0.025 - x[1] * x[1])', '0')
bcu_outflow = fe.DirichletBC(K.sub(1), fe.Expression(outflow_profile, degree = 2), outflow)
bcu_outflow_p = fe.DirichletBC(K.sub(0), fe.Constant(0), lower_boundary_pressure_fixed_point, method ='pointwise')
bcu = [bcu_inflow, bcu_wall, bcu_stress, bcu_wall_stress, bcu_outflow, bcu_outflow_p]
# Without surface pressure
F_st_NSE = rho * fe.dot(fe.dot(u, fe.nabla_grad(u)), v) * fe.dx \
- fe.div(v) * p * fe.dx \
+ fe.inner(stress(u, c), epsilon(v)) * fe.dx \
- fe.dot(fe.dot(stress(u,c), n), v) * fe.ds \
+ fe.div(u) * q * fe.dx \
# F_st_NSE = rho * fe.dot(fe.dot(u, fe.nabla_grad(u)), v) * fe.dx \
# - fe.div(v) * p * fe.dx \
# + fe.dot(p * n, v) * fe.ds \
# + fe.inner(stress(u, c), epsilon(v)) * fe.dx \
# - fe.dot(fe.dot(stress(u,c), n), v) * fe.ds \
# + fe.div(u) * q * fe.dx \
F_st_UCM = fe.inner(fe.dot(u, fe.nabla_grad(c)), d) * fe.dx \
- fe.inner(fe.dot(c, fe.nabla_grad(u)), d) * fe.dx \
- fe.inner(fe.dot(fe.grad(u), c), d) * fe.dx \
+ fe.inner((c-fe.Identity(2)), d) * fe.dx \
F = F_st_NSE + F_st_UCM
J_st = fe.derivative(F, z, z_)
problem = fe.NonlinearVariationalProblem(F, z, bcs=bcu, J=J_st)
solver = fe.NonlinearVariationalSolver(problem)
solver.solve()
p_, u_, c_ = z.split(deepcopy = True)
save_sol(c_, u_, p_, t)
u_prev.assign(u_)
c_prev.assign(c_)
for step in range(num_steps-1):
step += 1
# Update time
t += dt
# Print Current time step and time
if rank == 0:
print("Start time step{}, t={}".format(step+1, t))
sys.stdout.flush()
# Define boundary condition
inflow_profile = ('0.01', '0')
bcu_inflow = fe.DirichletBC(K.sub(1), fe.Expression(inflow_profile, degree =2), inflow)
bcu_wall = fe.DirichletBC(K.sub(1), fe.Constant((0,0)), walls)
bcu_stress = fe.DirichletBC(K.sub(2), fe.Constant([[0,0], [0,0]]), inflow)
bcu_wall_stress = fe.DirichletBC(K.sub(2), fe.Constant([[0,0], [0,0]]), walls)
outflow_profile = ('0.01 * 1600 * (0.025 * 0.025 - x[1] * x[1])', '0')
bcu_outflow = fe.DirichletBC(K.sub(1), fe.Expression(outflow_profile, degree = 2), outflow)
bcu_outflow_p = fe.DirichletBC(K.sub(0), fe.Constant(0), lower_boundary_pressure_fixed_point, method ='pointwise')
bcu = [bcu_inflow, bcu_wall, bcu_stress, bcu_wall_stress, bcu_outflow, bcu_outflow_p]
# Define finite derivative
du = der_be(u, u_prev, k)
dc = der_be_UCM(c, c_prev, k)
F_t = rho*fe.dot(du, v)*fe.dx \
+ rho * fe.dot(fe.dot(u, fe.nabla_grad(u)), v) * fe.dx \
- fe.div(v) * p * fe.dx \
- fe.inner(stress(u, c), epsilon(v)) * fe.dx \
+ fe.dot(fe.dot(stress(u,c), n), v) * fe.ds \
+ fe.div(u) * q * fe.dx \
F_t_UCM = fe.inner(dc, d) * fe.dx \
+ fe.inner(fe.dot(u, fe.nabla_grad(c)), d) * fe.dx \
- fe.inner(fe.dot(c, fe.nabla_grad(u)), d) * fe.dx \
- fe.inner(fe.dot(fe.grad(u), c), d) * fe.dx \
+ fe.inner((c-fe.Identity(2)), d) * fe.dx \
F_transient = F_t + F_t_UCM
J_t = fe.derivative(F_transient, z, z_)
problem = fe.NonlinearVariationalProblem(F_transient, z, bcs=bcu, J=J_t)
solver = fe.NonlinearVariationalSolver(problem)
solver.solve()
p_, u_, c_ = z.split(deepcopy = True)
save_sol(c_, u_, p_, t)
u_prev.assign(u_)
c_prev.assign(c_)
if rank == 0:
print("Elapsed time: ", time.time() - start)
Thank you!