Hi everyone,
I’m quite new to FEniCS and am simulating the Lid driven Cavity flow with viscoelastic term (UCM model). The mesh is 1 x 1 square.
Calculate the formulation with NonlinearvariationalSolver as it is shown below…
And I got error message that the newton residual value does not converge, diverges to -nan value.
What I found to solve this was that I need to give the initial value. And that’s why I gave the value of tensor, vector and pressure value.
Can you please help me to solve this problem?
Thank you for your answers.
Best regards,
Jo
import fenics as fe
import numpy as np
import os, sys
import matplotlib.pyplot as plt
from mpi4py import MPI
rank = MPI.COMM_WORLD.Get_rank()
tag = ‘LDC_viscoelastic’
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()
mesh_name = ‘mesh’
Import external mesh
mesh_dir = os.path.join(‘mesh’, mesh_name)
mesh = fe.Mesh(os.path.join(mesh_dir, mesh_name + ‘_mesh.xml’))
Define boundaries
boundary = fe.MeshFunction(‘size_t’, mesh, os.path.join(mesh_dir, mesh_name + ‘_facet.xml’))
Define Function spaces with mixed elements
P = fe.FiniteElement(‘CG’, fe.triangle, degree = 1)
V = fe.VectorElement(‘CG’, fe.triangle, degree = 2)
T = fe.TensorElement(‘CG’, fe.triangle, degree = 2)
elem = fe.MixedElement([T, V, P])
Q = fe.FunctionSpace(mesh, elem)
K = fe.VectorFunctionSpace(mesh, ‘CG’, 2)
O = fe.TensorFunctionSpace(mesh, ‘CG’, 2)
Define trial and test functions
z = fe.Function(Q)
z_ = fe.TrialFunction(Q)
c, u, p = fe.split(z)
d, v, q = fe.TestFunctions(Q)
Define expressions used in variational form
n = fe.FacetNormal(mesh)
Define a function to save u, p , stress and pressure gradient
def save_sol(c_, u_, p_, Reynolds):
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 Reynolds == 0:
flag = False
else:
flag = True
Save solutions to XDMF files
xdmf_u.write_checkpoint(u_, “u”, Reynolds, fe.XDMFFile.Encoding.HDF5, flag)
xdmf_p.write_checkpoint(p_, “p”, Reynolds, fe.XDMFFile.Encoding.HDF5, flag)
xdmf_s.write_checkpoint(strs, “s”, Reynolds, fe.XDMFFile.Encoding.HDF5, flag)
xdmf_dp.write_checkpoint(dp, “dp”, Reynolds, fe.XDMFFile.Encoding.HDF5, flag)
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 variational Problem
Define symmetric gradient, sigma(u,p) and stress(u)
def epsilon(u):
return fe.sym(fe.nabla_grad(u))
def stress(u, c):
return G * (c - fe.Identity(len(u)))
Define boundary condition
top_velocity = 1
bc_inflow = fe.DirichletBC(Q.sub(1), fe.Constant((1, 0)), boundary, 1)
bc_noflow1 = fe.DirichletBC(Q.sub(1), fe.Constant((0, 0)), boundary, 2)
bc_noflow2 = fe.DirichletBC(Q.sub(1), fe.Constant((0, 0)), boundary, 3)
bc_noflow3 = fe.DirichletBC(Q.sub(1), fe.Constant((0, 0)), boundary, 4)
bc_pressure = fe.DirichletBC(Q.sub(2), fe.Constant(0), “near(x[0], -0.5) && x[1] == -0.5”,“pointwise”)
boundaries = [bc_inflow, bc_noflow1, bc_noflow2, bc_noflow3, bc_pressure]
Constants
G = 1
mu = 0.1
rho = 1
G = fe.Constant(G)
mu = fe.Constant(mu)
rho = fe.Constant(rho)
Reynolds = rho * top_velocity / mu
Reynolds = fe.Constant(Reynolds)
print(“Doing it for Re = {}”.format(Reynolds))
Initial Value
e_c0 = fe.Expression(((‘0.0’, ‘0.0’), (‘0.0’,‘0.0’)), degree = 2)
e_u0 = fe.Expression((‘0.0’, ‘0.0’), degree = 2)
e_p0 = fe.Expression(‘0.0’, degree = 1)
p0 = fe.interpolate(e_p0, Q.sub(2).collapse())
u0 = fe.interpolate(e_u0, Q.sub(1).collapse())
c0 = fe.interpolate(e_c0, Q.sub(0).collapse())
fe.assign(z, [c0, u0, p0])
Variational Formulation
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
+ G * fe.inner(stress(u, c), epsilon(v)) * fe.dx
- G * 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(len(u))), d) * fe.dx
F_st = F_st_NSE
+ F_st_UCM
Build Jacobian Matrix
J_st = fe.derivative(F_st, z, z_)
Define NonlinearVaritationalProblem and solve
problem = fe.NonlinearVariationalProblem(F_st, z, bcs = boundaries, J = J_st)
solver = fe.NonlinearVariationalSolver(problem)
solver.solve()
c_, u_, p_ = z.split(deepcopy = True)
Save solution
save_sol(c_, u_, p_, Reynolds)
The error message
Blockquote
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.017e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 2: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 3: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f