Whether the solver Settings are reasonable

General
1

Hi everyone, my question has already been addressed here:https://fenicsproject.discourse.group/t/the-solution-speed-is-very-slow/11576/3. I have modified my code according to your suggestion as follows:

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, Can you tell me if it is normal:

1687693064364
1687693064364727×506 12.7 KB

1687693533168(1)

2

But the solution will stop automatically after a certain time(win+docker+fenics2019).

3

Its normal output from solving a non-linear problem.
What do you mean when you say:

Do you Get an error message?

4

The calculation stops without receiving any error message, and the notebook loses connection

5

Try to run it as a Python script and not a notebook, as it usually captures more errors.

6

I run the following command on the terminal to solve for the.py file:
docker run --rm -v E:\fenics_work:/home/fenics/shared -w /home/fenics/shared quay.io/fenicsproject/stable "python3 heat_diffsion.py"
However, when the solution reaches 70%, it always stops the calculation, but does not suffer any errors, which seems to be caused by docker crashes.

7

Please provide a minimal (complete) reproducible code; that means a code that can be copy-pasted by me or anyone else into an editor, and be executed directly without modifications.

I would also strongly recommend to use the image:
ghcr.io/scientificcomputing/fenics:2023-04-21
as it is built on an up to date linux machine.

8

Here’s 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)
9

My docker running memory is 4GB, when the solving progress reaches 30%, the occupied memory is 3.8G, and then when the solving progress reaches 70%, the calculation will stop without any error, and docker will also be stuck and unable to refresh.

10

Your code has inefficiencies. See, e.g. Out of memory error in frequency for loop - #4 by nate.

1 Like
11

Thank you, this is really the problem here (this problem you have proposed to me), but I am a beginner, do not know how to set up the iterative solver correctly, can you recommend some materials for me to learn the relevant theory?

12

I already showed you this in: The solution speed is very slow - #3 by ss780060872

Which avoids recreation of matrices, and recompilation of forms.