Converting 3D heat equation to a x - y plot with respects to heat and time

Hello,

I am trying to make a simple x - y graph to represent heat vs time with my cylindrical heat equation. I am relatively new to fenics/python and would please like some help.

Thanks,

My code and Results:

from future import print_function
from dolfin import *
from mshr import *
from fenics import *
import time

#Making a cylindrical geometry
cylinder = Cylinder(Point(0, 0, -1), Point(0, 0, 1), 0.5, 0.5)
geometry = cylinder

Making Mesh (40 corresponds to the mesh density)

mesh = generate_mesh(geometry, 40)

Save the mesh

File(“cylinder.pvd”) << mesh

V = FunctionSpace(mesh, ‘P’, 1)

T = 2.0 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size

Define boundary condition

def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

Define initial value

u_0 = Expression(‘exp(-apow(x[0], 2) - apow(x[1], 2))’,
degree=2, a=5)
u_n = interpolate(u_0, V)

Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

F = uvdx + dt*dot(grad(u), grad(v))dx - (u_n + dtf)vdx
a, L = lhs(F), rhs(F)

Create VTK file for saving solution

vtkfile = File(‘heat_gaussian/solution.pvd’)

Time-stepping

u = Function(V)
t = 0
for n in range(num_steps):

# Update current time
t += dt

# Compute solution
solve(a == L, u, bc)

# Save to file and plot solution
vtkfile << (u, t)
plot(u)

# Update previous solution
u_n.assign(u)

import matplotlib.pyplot as plt
plt.show()

You can use paraview for this,which has the option of creating slices.

1 Like

You can make an empty list for heat and time values. For example put time_vals = [] before your time-stepping loop and then store the time values at each time step in your loop by using append time_vals.append(t)

Repeat this process to store heat values with a new list heat_vals = []. You will then have time and heat lists which you can scatterplot with import matplotlib.pyplot as plt

plt.xlabel('Time')
plt.ylabel('Heat')
plt.ylabel('Heat vs Time')
plt.scatter(time_vals,heat_vals)
plt.show()
1 Like

Thanks, where should I place the heat_vals = []? and should it’s append also contain t?

Thanks, I tried running it in paraview but I am getting the error below. I downloaded the pdv file onto my desktop and tried applying it but got the error in the screenshot. Any ideas?

Define the empty list outside the loop. You want to append heat values to the list at each time step in the loop. You can just append u(t) to the list at each step in the loop, since u(t) is your heat at time t I’m presuming?

1 Like

I would try paraview 5.7, as there was a bug in 5.8.0, and 5.8.1 is not an official release yet

3 Likes

I second this. 5.7 is stable, 5.8 is broken for now.

1 Like