Hi everybody,

I’ve recently started working with Fenics and I have a problem with the final solution which is a distribution of temperature in the steady-state case; That is, when we don’t have any temporal variation and the system has reached an equilibrium. The problem is that I notice a wild oscillation throughout the border where I forced the temperature to remain 20 degrees C by using *Outer_bc = DirichletBC(V, 20, DomainBoundary())*. I did a simple simulation for a cylinder where I expect to see a pretty uniform distribution of 37 degrees C in the middle and a gradual decrease to 20 C towards the border which has a fixed temperature (please see the snippet of my code for this scenario).

First I had used a 1st order CG for the function space, so I suspected the problem might have come from poor meshing. So I changed the function space to 2nd order CG and it seemed to solve the problem (please see the attached temperature profiles of the two simulations). I was happy and I thought the problem would completely disappear as I kept increasing the resolution by importing a finer and finer mesh. So I decided to make my mesh even finer to get a smoother profile and I didn’t expect to see the old problem resurface. But to my consternation, it did. I increased the number of mesh elements from 400K to almost 4 million and sure I ran out of memory when I tried to run the code in series. So I ran it in parallel and in the final result I got very good and smooth distribution in the middle but those overshoots to 55 degrees C and even higher came back around the edges!

Is this a problem of meshing really or am I missing something here? To me, it doesn’t seem to come from just poor meshing because it happened even for a very fine mesh at the end, but why didn’t it happen for the 2nd run when I just used the old (coarser) mesh with just increased degrees of freedom (2nd order CG)??

I’m wondering if someone has experienced the same problem. It’s highly appreciated if you share your thought and experiences in this regard and many thanks in advance.

`from fenics import * #%% Declaration of Constants Kapa_value = [0.5, 0.55] #Thermal Conductivity (W/m/'c), [Muscle, Tumor] Cp_value = [3421, 3500]#Heat Capacity (J/Kg/'c), [Muscle, Tumor] Cb = 3617 #specific heat (Heat Capacity) of blood (J/Kg/'c) Wb_value = [0.8, 0.5] #local perfusion rate (kg/m3/s), [Muscle, Tumor] Tb = 37 #Arterial blood temperature Rho_value = [1090, 1040]#Density of Muscle (Kg/m3) Bolus_Temp = Constant(20) Body_Temp = Constant(37) #Normal body temperature #%% Classic Mesh conversion using dolfin-convert (legacy)! mesh = Mesh('../GeoMesh.xml') Markers = MeshFunction('size_t',mesh,'../GeoMesh_physical_region.xml') BCs = MeshFunction('size_t',mesh,'../GeoMesh_facet_region.xml') dx = Measure('dx', domain=mesh, subdomain_data=Markers) #%% Subdomains' property assignment V0 = FunctionSpace(mesh, 'DG', 0) Kapa = Function(V0) Cp = Function(V0) Wb = Function(V0) Rho = Function(V0) for n, cell_no in enumerate(Markers.array()): Subdomain_Tag = int(cell_no - 1) #Markers.array()[n] Kapa.vector()[n] = Kapa_value[Subdomain_Tag] Cp.vector()[n] = Cp_value[Subdomain_Tag] Wb.vector()[n] = Wb_value[Subdomain_Tag] Rho.vector()[n] = Rho_value[Subdomain_Tag] MeshExport = File('CylinderWithoutF/mesh.pvd') MeshExport << mesh ThermalCond = File('CylinderWithoutF/Kapa.pvd') ThermalCond << Kapa #%% Boundary Condition V = FunctionSpace(mesh, 'CG', 2) Outer_bc = DirichletBC(V, Bolus_Temp, DomainBoundary()) #%% Define Variational Problem T = TrialFunction(V) w = TestFunction(V) #%% Symbolic Bioheat equation F = -1*Kapa*dot(grad(T), grad(w))*dx + (dot(grad(T), grad(Kapa)) - Cb*Wb*(T - Tb))*w*dx a, L = lhs(F), rhs(F) # Create VTK file for saving solution vtkfile = File('CylinderWithoutF/Temperature.pvd') # Export back to Matlab hdfHandle = HDF5File(mesh.mpi_comm(), 'CylinderWithoutF/FinalTemperature.h5','w') hdfHandle.write(mesh, '/mesh') # Compute the solution T = Function(V) #solve(a == L, T, Outer_bc) solve(a == L, T, Outer_bc, solver_parameters = {'linear_solver':'gmres'}) # Save to file and plot solution vtkfile << T hdfHandle.write(T, '/Temperature')`

Here are the results visualized in Paraview [from left to right: CG1 for 400K elements, CG2 for 400K, and finally CG2 for 4Million mesh elements!) Please notice the high red zones stretched across the border which are the problematic parts and numerical artifacts and they shouldn’t be there.