# The solution speed is very slow

Hello everyone, I encountered the following three problems when solving nonlinear heat transfer problems(windows+docker+fenics2019.1.0):
1.When I change the total solving time, it always doesn’t converge, I usually refine the grid and reduce the time step to make it converge, but this tends to increase the time cost, how do I optimize my code?
2.the resulting.h5 file is very large, about 7GB,Is that normal?
3.How can I detect iterations or progress more closely?
here is my code:

``````from fenics import *
import matplotlib.pyplot as plt
import numpy as np

L = 0.01
W = 0.01
nx = 200
ny = 200
mesh = RectangleMesh(Point(0, 0), Point(L, W), nx, ny)
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=0)
#h = Expression('70000', degree=0)
T_ini = Expression('298.15', degree=0)
#T_env = Expression('323.15', degree=0)
f = Constant(0.0)

t_total = 100
num_steps = 5000
dt = t_total / num_steps

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

class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x-0.01)<DOLFIN_EPS

class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x-0.00)<DOLFIN_EPS

class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x-0.01)<DOLFIN_EPS

bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()

bc_bottom = DirichletBC(V, T_ini, by0)

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

#boundary_markers.set_all(9999)
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)

boundary_conditions = {0: {'Neumann': Constant(0)}, 1: {'Neumann': Constant(0)}, 2: {'Dirichlet': T_ini}, 3: {'Neumann': q}}

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]:
#if boundary_conditions[i]['Neumann'] != 0:
q = boundary_conditions[i]['Neumann']
integrals_N.append(q*v*ds(i))

integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
h, T_envi = boundary_conditions[i]['Robin']
integrals_R.append(h*(T - T_envi)*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)

xdmffile_T = XDMFFile('10_10_100s_slab/temperature.xdmf')

timeseries_T = TimeSeries('10_10_100s_slab/temperature_series')

File('10_10_100s_slab/mesh.xml.gz') << mesh

xdmffile_T.parameters["flush_output"] = True
xdmffile_T.parameters["functions_share_mesh"] = True

progress = Progress('Time-stepping', num_steps)

t = 0
times = 
temperatures=[298.15]
plt.figure(figsize=(12,8))
for n in range(num_steps):

t += dt

solve(F == 0, T, bcs)
plot(T)

point = Point(0.005, 0.01)
temperature = T(point)
temperatures.append(temperature)
times.append(t)

xdmffile_T.write(T, t)

timeseries_T.store(T.vector(), t)

T_n.assign(T)

set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
``````

Be grateful!

Avoid calling solve F==0 inside the loop, as it recreate a lot of structures (sparse matrices, vector for thr RHS, ksp solver etc.

This is covered in various other posts on the forum. In the context of DOLFINx at:

and for legacy dolfin and linear problems in

You are not customizing your solver. See for instance:

You can write a custom newton solver, see for instance Set krylov linear solver paramters in newton solver - #3 by nate

Hello dokken, thank you very much for your answer. I have modified my code according to your suggestion as follows:

``````progress = Progress('Time-stepping', num_steps)
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-10
prm['newton_solver']['relative_tolerance'] = 1e-9
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['linear_solver'] = 'mumps'

t = 0
plt.figure(figsize=(12,8))
for n in range(num_steps):

t += dt

solver.solve()
plot(T)

xdmffile_T.write(T, t)

timeseries_T.store(T.vector(), t)

T_n.assign(T)

set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
set_log_level(LogLevel.INFO)
``````

I am a beginner, I may not know the specific meaning of the above solver, but I can solve normally, and the speed is much faster, the following is some information output by the terminal, please help me to see if it is normal: You are still using a direct solver for the linearized problem, which is not as efficient as using an iterative one.

Maybe my theory on this aspect is wrong, please forgive my stupidity, could you give me some specific guidance?

I would suggest looking at: How to choose the optimal solver for a PDE problem? - #2 by nate
and
Default absolute tolerance and relative tolerance - #4 by nate
as it gives a lot of general guidance.

If you want to understand what happens in a “Newton solve” i would suggest reading Custom Newton solvers — FEniCSx tutorial