First, it seems like you have a scaling mistake in your code:
What happened to dt for the ds term.
Secondly, your code can be simplified quite alot:
mesh, cell_markers, facet_markers = gmshio.read_from_msh("the_mesh.msh", MPI.COMM_WORLD, gdim=3)
# Define function space
V = FunctionSpace(mesh, ("Lagrange", 2))
# Define the "unknown" variable of the equation.
Temp = Function(V)
# Physical parameters.
κ = 49.0
C_p = 120
density = 9.7e3 # kg / m^3.
T_amb = 0.0
q_SB = 5.6704e-8
h = 10
T_fixed_surface = 21
convection_surface = 22
T_left = 100
# Define Dirichlet boundary conditions to the facets corresponding to the left end of the rod.
left_facets = facet_markers.find(T_fixed_surface)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bc = dirichletbc(ScalarType(T_left), left_dofs, V)
bcs = [bc]
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh,subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
# Initial temperature of the rod.
T_0 = 0.0
def temp_init(x):
values = np.full(x.shape[1], T_0, dtype = ScalarType)
return values
Temp.name = "Temp"
Temp.interpolate(temp_init)
# Simulation parameters.
t = 0 # Start time
t_final = 6.0 # Final time
num_steps = 100
dt = t_final / num_steps # time step size
T_n = Function(V)
T_n.name = "T_n"
T_n.interpolate(temp_init)
v = TestFunction(V)
F_T = density * C_p * (Temp - T_n )* v * dx + dt * dot(κ * grad(Temp), grad(v)) * dx + dt * v * h * (Temp - T_amb) * ds(convection_surface)
# Define the problem.
problem = NonlinearProblem(F_T, Temp, bcs=bcs)
solver = NewtonSolver(mesh.comm, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-14
solver.report = True
ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000
ksp.setFromOptions()
times, Qs_ratios, spatial_average_temperatures = [], [], []
n = FacetNormal(mesh)
left_flux = form(dot(-κ * grad(Temp), n)*ds(T_fixed_surface))
right_flux = form(dot(-κ * grad(Temp), n)*ds(convection_surface))
Q = VectorFunctionSpace(mesh, ("DG", 1))
q = Function(Q)
flux_calculator = Expression(-κ * grad(Temp), Q.element.interpolation_points())
out_file = VTXWriter(mesh.comm, "t.bp", [Temp], engine="BP4")
out_flux = VTXWriter(mesh.comm, "flux.bp", [q], engine="BP4")
#set_log_level(LogLevel.INFO)
for time_step in np.linspace(t, t_final, num_steps):
times.append(t)
# Compute total heat convected away.
# Heat entering the rod from the left side.
Q_in = assemble_scalar(left_flux)
# Heat leaving the rod on all other surfaces.
Q_out = assemble_scalar(right_flux)
if Q_in == 0:
Qs_ratio = 0
else:
Qs_ratio = abs(Q_out) / abs(Q_in)
print(Q_out, Q_in, time_step)
Qs_ratios.append(Qs_ratio * 100)
t += dt
out_file.write(t)
n, converged = solver.solve(Temp)
T_n.x.array[:] = Temp.x.array
q.interpolate(flux_calculator)
out_flux.write(time_step)
#spatial_average_temperature = assemble_scalar(form(Temp * dx(domain = mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))
#spatial_average_temperatures.append(spatial_average_temperature)
out_flux.close()
out_file.close()
#plt.plot(times, spatial_average_temperatures)
plt.plot(times, Qs_ratios)
plt.savefig("Test.png")
plt.close()
I would also like to note that your mesh is incredibly coarse (and it is scaled such that the cell jacobian is very tiny, which could lead to numerical issues. For instance the fluxes are very large (1e6 magnitude)
You are also not attaching the appropriate null space (the constant space) to your multi-grid solver.
Considering the code above with P1, P2, P3 gives
