NAN occurs when a custom grid is computed

General
1

Hi everyone, I currently create a 1D grid using numpy generated data, but when solving a heat transfer problem prompts NAN, the code is as follows:

vertices = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])

mesh = Mesh()
editor = MeshEditor()

editor.open(mesh, 'interval', 1, 1)
editor.init_vertices(len(vertices))  
editor.init_cells(len(vertices)-1)    

for i in range(len(vertices)):
    editor.add_vertex(i, np.array([vertices[i]]))
    
for j in range(len(vertices)-1):
    editor.add_cell(j, np.array([j, j+1]))

editor.close()

V = FunctionSpace(mesh, 'P', 1)

.............

class BoundaryX0(SubDomain):  
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.00) < DOLFIN_EPS

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 1) < DOLFIN_EPS

bx0 = BoundaryX0()  
bx1 = BoundaryX1()

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)  

bx0.mark(boundary_markers, 0)   
bx1.mark(boundary_markers, 1)

boundary_conditions = {0: {'Neumann': q}, 1: {'Dirichlet': T_ini}}

bcs = []   
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
        bcs.append(bc)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
T_n = interpolate(T_ini, V)

T = Function(V)
v = TestFunction(V)

integrals_N = []    
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        q = boundary_conditions[i]['Neumann']
        integrals_N.append(q*v*ds(i))

F = mu*c*T*v*dx + dt*k(T)*dot(grad(T), grad(v))*dx - mu*c*T_n*v*dx - dt*f*v*dx - dt*sum(integrals_N)

.............

Thank you very much!

2

Without a working example it’s difficult to give advice. For example we don’t know what k(T), mu and c are.

3

Here is the complete code:

vertices = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])

mesh = Mesh()
editor = MeshEditor()

editor.open(mesh, 'interval', 1, 1)
editor.init_vertices(len(vertices))  
editor.init_cells(len(vertices)-1)    

for i in range(len(vertices)):
    editor.add_vertex(i, np.array([vertices[i]]))
    
for j in range(len(vertices)-1):
    editor.add_cell(j, np.array([j, j+1]))

editor.close()

V = FunctionSpace(mesh, 'P', 1)

def k(T):  
    return 174.9274 - 0.1067 * T + 5.0067e-5 * T**2 - 7.8349e-9 * T**3

mu = Constant(19300)  
c = Constant(137)    
q = Expression('1e7', degree=1)   
T_ini = Expression('298.15', degree=1)    
f = Constant(0.0)  

t_total = 10   
num_steps = 500    
dt = t_total / num_steps  

class BoundaryX0(SubDomain):   
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.00) < DOLFIN_EPS

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 1) < DOLFIN_EPS

bx0 = BoundaryX0() 
bx1 = BoundaryX1()

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1) 

bx0.mark(boundary_markers, 0) 
bx1.mark(boundary_markers, 1)

boundary_conditions = {0: {'Neumann': q}, 1: {'Dirichlet': T_ini}}

bcs = []  
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
        bcs.append(bc)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
        
T_n = interpolate(T_ini, V)

T = Function(V)
v = TestFunction(V)

integrals_N = []  
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        q = boundary_conditions[i]['Neumann']
        integrals_N.append(q*v*ds(i))

F = mu*c*T*v*dx + dt*k(T)*dot(grad(T), grad(v))*dx - mu*c*T_n*v*dx - dt*f*v*dx - dt*sum(integrals_N)

J = derivative(F, T)
problem = NonlinearVariationalProblem(F, T, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['absolute_tolerance'] = 1e-6
prm['newton_solver']['relative_tolerance'] = 1e-7
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0

t = 0
plt.figure(figsize=(12,8))
for n in range(num_steps):
    
    t += dt    
  
    solver.solve()
      
    T_n.assign(T)
4

Looks like you have a duplicate point at x=0.5 in

vertices = np.concatenate([np.linspace(0, 0.5, num=100),np.linspace(0.5, 1, num=200)])

which gives you a degenerate cell with no volume. If I account for this the system is solvable.