Can someone please help me? I don’t understand why the solution diverges, if you put viscosity_0 = const, the program works.
My code:
%%
from fenics import *
from mshr import *
import time
import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
%%
comm = MPI.COMM_WORLD
print(“Hello from process”, comm.rank)
%%
xdmffile_u = XDMFFile(comm,‘final_v_e_3/velocity.xdmf’)
xdmffile_p = XDMFFile(comm,‘final_v_e_3/pressure.xdmf’)
xdmffile_e = XDMFFile(comm,‘final_v_e_3/energy.xdmf’)
xdmffile_nu = XDMFFile(comm,‘final_v_e_3/viscos.xdmf’)
%%
%%
Create mesh
rect_width = 0.0018
rect_height = 0.016
trap_bottom_width = 0.0018
trap_top_width = 0.0004
trap_height = 0.0015
rectangle = Rectangle(Point(0, 0), Point(rect_width, rect_height))
trapezoid_bottom_left = Point(0, rect_height)
trapezoid_bottom_right = Point(trap_bottom_width, rect_height)
trapezoid_top_left = Point((trap_bottom_width - trap_top_width) / 2, rect_height + trap_height)
trapezoid_top_right = Point((trap_bottom_width + trap_top_width) / 2, rect_height + trap_height)
trapezoid = Polygon([trapezoid_bottom_left, trapezoid_bottom_right, trapezoid_top_right, trapezoid_top_left])
additional_rect_width = trap_top_width
additional_rect_height = 0.0005
additional_rectangle = Rectangle(Point((trap_bottom_width - trap_top_width) / 2, rect_height + trap_height),
Point((trap_bottom_width + trap_top_width) / 2, rect_height + trap_height + additional_rect_height))
domain = rectangle + trapezoid + additional_rectangle
mesh = generate_mesh(domain, 128)
num_steps = 1000 # number of time steps
dt = 1/num_steps # time step size
rho = 1122.80
rho = Constant(rho)
f = Constant((0, 0))
k = Constant(dt)
C_p = Constant(2140)
k_p = Constant(0.180)
n_n = 0.3846
tau = 129.0
D1 = 2.05e7
A1 = 16.71
T_sharp = 373.15
A2 = 51.60
theta=0.5
V = VectorElement(“Lagrange”,triangle, 2)
Q = FiniteElement(“Lagrange”,triangle, 1)
E = FiniteElement(“Lagrange”,triangle, 1)
W_element = MixedElement(V, Q,E)
W = FunctionSpace(mesh, W_element)
NU = FunctionSpace(mesh, ‘P’, 1)
nu = Function(NU)
y = TestFunction(W)
(v,q,g) = split(y)
w0 = Function(W)
(u0,p0,e0) = split(w0)
w = Function(W)
(u,p,e) = split(w)
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0)
class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[1], rect_height + trap_height + additional_rect_height))
class Walls(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and not (near(x[1], 0) or (near(x[1], rect_height + trap_height + additional_rect_height)))
boundaries = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1, 0)
inlet = Inlet()
outlet = Outlet()
walls = Walls()
inlet.mark(boundaries, 1)
outlet.mark(boundaries, 2)
walls.mark(boundaries, 3)
inflow = Expression((‘0’, ‘0.005’,), degree=0)
noslip = Constant((0.0, 0.0))
Temperature = Constant(468)
bc_inlet_v = DirichletBC(W.sub(0), inflow, boundaries, 1)
bc_walls = DirichletBC(W.sub(0), noslip, boundaries, 3)
bc_outlet = DirichletBC(W.sub(1), Constant(0), boundaries, 2)
bce_inlet = DirichletBC(W.sub(2), Constant(T_sharp + 2.0), boundaries, 1)
bce_walls = DirichletBC(W.sub(2),Temperature, boundaries, 3)
bce = [bce_inlet, bce_walls]
bc_v_p = [bc_inlet_v, bc_walls, bc_outlet]
bcs = bc_v_p + bce
def epsilon(u):
return sym(nabla_grad(u))
def viscosity_0(e_):
return D1exp(-A1(e_ - Constant(T_sharp))/(Constant(A2)+(e_ - Constant(T_sharp))))
def viscosity(w, e):
gamma_05 = pow(0.5*tr(epsilon(w)*epsilon(w)) + DOLFIN_EPS, 0.5)
return viscosity_0(e)/(1+ pow((viscosity_0(e)*gamma_05)/tau, 1. - n_n))
def c(v,e,g) :
return (dot(k_pgrad(e),grad(g)) + C_prho*dot(v,grad(e))g - 2viscosity(v,e)*tr(epsilon(v)*epsilon(v))*g)*dx
F_v_p = rho*dot(dot(u, nabla_grad(u)), v)dx
+ 2viscosity(u,e)*inner(epsilon(u), nabla_grad(v))*dx
- div(v)pdx - div(u)qdx
F_e = c(u, e, g)
F_e_0 = c(u0, e0, g)
F_v_p_0 = rho*dot(dot(u0, nabla_grad(u0)), v)dx
+ 2viscosity(u0,e0)*inner(epsilon(u0), nabla_grad(v))*dx
- div(v)p0dx - div(u0)qdx
F = rhoinner((u - u0) / k, v)dx + C_prho(inner((e-e0),g)/k)dx +(1.0-theta)(F_v_p_0 +F_e_0)+ theta*(F_v_p + F_e)
J = derivative(F, w)
problem = NonlinearVariationalProblem(F,w,bcs,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm[‘nonlinear_solver’] = ‘newton’
prm[‘newton_solver’][‘absolute_tolerance’] = 1E-12
prm[‘newton_solver’][‘relative_tolerance’] = 1e-12
t = dt
start = time.time()
for n in range(5000):
begin(“Solving …”)
solver.solve()
end()
xdmffile_u.write(w.sub(0), t)
xdmffile_p.write(w.sub(1), t)
xdmffile_e.write(w.sub(2), t)
nu.assign(project(viscosity(w.sub(0), w.sub(2)), NU))
xdmffile_nu.write(nu, t)
assign(w0,w)
end = time.time()
print(“CPU(s)”, end - start)
error:
Process 0: Solving …
Process 0: Solving nonlinear variational problem.
Process 1: Solving …
Process 1: Solving nonlinear variational problem.
Process 3: Solving …
Process 3: Solving nonlinear variational problem.
Process 2: Solving …
Process 2: Solving nonlinear variational problem.
Process 0: Newton iteration 0: r (abs) = 3.950e+04 (tol = 1.000e-12) r (rel) = 1.000e+00 (tol = 1.000e-12)
Process 0: Newton iteration 1: r (abs) = nan (tol = 1.000e-12) r (rel) = nan (tol = 1.000e-12)
Traceback (most recent call last):
File “/home/artlinux1998/project/rd_3.py”, line 218, in
solver.solve()
RuntimeError:
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 2