The solution speed is very slow

[ss780060872]

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]-0.00)<DOLFIN_EPS

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

class BoundaryY1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1]-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 = [0]
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!

[dokken]

@ss780060872

for n in range(num_steps):
    
    t += dt
    
    solve(F == 0, T, bcs)

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:

Out of memory error in frequency for loop

dolfinx.fem.Constant wraps constant quantities which can change between assemblies of the underlying matrix. dolfinx.fem.petsc.LinearProblem has a PETSc Mat and Vec. Creating one every time step is extremely expensive, and not their intended use. ... # Problem setup k = dolfinx.fem.Constant(mesh, 0.0) k.value = some_important_value ... # Weak form definition problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) ... # Other code for nf in range(0,len(f_a…

and for legacy dolfin and linear problems in

Simula-SSCP/SSCP_2023_lectures/blob/main/L12 (FEniCS Intro)/L02_FEniCS_Diffusion.ipynb

@ss780060872

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?

You are not customizing your solver. See for instance:

Default absolute tolerance and relative tolerance

I write the following in the context of FE modelling rather than numerical analysis. I.e. this should be interpreted as tips and tricks. Perhaps someone else can contribute something more rigorous, or their own experience. Do you know what to expect the residual of your system to be to machine precision? Adjust your absolute tolerance accordingly. E.g: A Poisson problem with a diffusion coefficient of 1 may yield a residual to machine precision around 10^{-12}. Linear elasticity with a Young’…

@ss780060872

3.How can I detect iterations or progress more closely?
here is my code:

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

[ss780060872]

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:

[dokken]

You are still using a direct solver for the linearized problem, which is not as efficient as using an iterative one.

[ss780060872]

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

[dokken]

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