I am trying to solve the non-linear transient heat equation.
When I execute the below MWE, the resultant temperature field is discontinuous (see picture attached). However, when the thermal properties (rho, Cv, conductivity) are set equal to 1, the resultant temperature field is continuous (see picture attached).
Could anyone please offer any advice as to why the PDE becomes ill-conditioned with these parameters? and how I could fix it?
from fenics import * import mshr #Define unstructured example mesh domain = mshr.Box(Point(0,0,0), Point(1,1,1)) mesh = mshr.generate_mesh(domain, 40) #Boundary condition for heat flux class S1(SubDomain): def inside(self, x, on_boundary): return on_boundary and near(x,0) S1 = S1() Surfaces = MeshFunction('size_t', mesh, 2) Surfaces.set_all(0) S1.mark(Surfaces, 1) ds = Measure("ds")(subdomain_data=Surfaces) #Define X-functions A = FunctionSpace(mesh, "CG", 1) Tn = interpolate(Constant(0), A) Tn_plus_1 = Function(A) vT = TestFunction(A) duT = TrialFunctions(A) #Thermal parameters rho = 3000 Cv = 800 conductivity = 20 timeStep = 1 for i in range(2): #Weak form of heat equation thermalForm = rho*Cv*Tn_plus_1*vT*dx - rho*Cv*Tn*vT*dx + Constant(timeStep)*inner(conductivity*grad(Tn_plus_1),grad(vT))*dx - Constant(timeStep)*inner(Constant(1e6),vT)*ds(1) #Construct non-linear variational problem J = derivative(thermalForm, Tn_plus_1, duT) problem = NonlinearVariationalProblem(thermalForm, Tn_plus_1, J=J) solver = NonlinearVariationalSolver(problem) solver.solve() #Assign current solution Tn.assign(Tn_plus_1) #Save time step to file File('Visualisation/T' + str(i) + '.pvd') << Tn_plus_1