Solving the weak formulation of the heat problem for a hexahedron cube

General
1

Heat Problem Weak Form Source of Hexahedron Cube:


import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
import numpy as np
import basix
from petsc4py import PETSc
from dolfinx.fem import Function, locate_dofs_topological, apply_lifting, set_bc
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.io import XDMFFile

# Define the material properties
κ = 79.5  # Thermal conductivity (W/(m·K))
ρ = 7850  # Density of the material (kg/m^3)
c_p = 0.451  # Specific heat capacity (J/(kg*K))
f = 50  # Heat source term (W/m^3)

# Create a unit cube mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 4, 4, 4, cell_type=dolfinx.mesh.CellType.hexahedron)

# Define the function space using your method
element_family = basix.ElementFamily.P
element_degree = 1
variant = basix.LagrangeVariant.equispaced
ufl_element = basix.ufl.element(element_family, mesh.topology.cell_name(), element_degree, variant)
V = dolfinx.fem.functionspace(mesh, ufl_element)

# Define the boundary condition function
u_D = dolfinx.fem.Function(V)
u_D.interpolate(lambda x: np.zeros(x.shape[1]))

# Locate boundary facets
def boundary(x):
    return np.logical_or.reduce((
        np.isclose(x[0], 0), np.isclose(x[0], 4),
        np.isclose(x[1], 0), np.isclose(x[1], 4),
        np.isclose(x[2], 0), np.isclose(x[2], 4)
    ))

boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, boundary)

# Locate the degrees of freedom
dofs = locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)

# Create a DirichletBC object
bc = dolfinx.fem.dirichletbc(u_D, dofs)
dirichlet_bc = dolfinx.fem.DirichletBC(bc)


# Define the test and trial functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define the time step and the initial condition
dt = 0.01  # Time step
u_n = Function(V)  # Initial condition (zero by default)
u_n.interpolate(lambda x: np.full(x.shape[1], 0))


# Define the variational problem
F = (ρ * c_p * u / dt) * v * ufl.dx - (ρ * c_p * u_n / dt) * v * ufl.dx + κ * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - f * v * ufl.dx

# Define the bilinear and linear forms
a = ufl.lhs(F)
L = ufl.rhs(F)

# Assemble the system matrix
A = assemble_matrix(dolfinx.fem.form(a))
A.assemble()
u = Function(V)

# Time-stepping
t = 0
T = 1  # Final time

def boundary_condition(x, t):
    return t * np.ones(x.shape[1])

# Create a XDMF file for output
with XDMFFile(MPI.COMM_WORLD, "output.xdmf", "w") as outfile:
    # Time-stepping
    t = 0
    T = 1  # Final time

    while t < T:
        u_D.interpolate(lambda x: boundary_condition(x, t))
        # Assemble the right-hand side
        b = assemble_vector(dolfinx.fem.form(L))

        # Apply boundary conditions to the vector
        apply_lifting(b, [dolfinx.fem.form(a)], [[bc]])
        set_bc(b, [bc])

        # Solve the system
        solver = PETSc.KSP().create(mesh.comm)
        solver.setOperators(A)
        solver.setType(PETSc.KSP.Type.PREONLY)
        solver.getPC().setType(PETSc.PC.Type.LU)
        solver.solve(b, u.vector)    

        # Update the solution for the next time step
        u_n.x.array[:] = u.x.array
        print(u.x.array)
        t += dt   

        # Write the solution to file
        outfile.write_function(u, t)
print('done...')

I did the best I could and have some type of solution coming out of the system. The units of the entry parameters are notated at the top. I hope those are OK…

Solution @ t=0:


      -8.57579853e-05,  1.77271103e-04, -8.57579853e-05,  1.77271103e-04,
       -1.04543838e-04,  2.15939191e-04, -8.57579853e-05,  1.77271103e-04,
       -9.98261863e-05,  2.06239969e-04,  1.09341942e-04, -9.98261863e-05,

Solution @ t=1:


[  3.26591142  -6.71856831  -6.71856831  13.85271957  -6.71856831
  13.85271957  13.85271957 -28.58819086  -3.9929514    8.24254544

So far the solutions I am seeing look a bit odd to me in terms of realistic values in terms of the negative values I am seeing so far. I don’t quite know yet why negation of temperatures might be happening. Does that seem in the norm?

So far I made some attempt at setting Dirichlet for such a problem but don’t really feel sure if that was done correctly.

So I would expect the solution here of this type of problem to be a temperature T of some sort. Sorry… little bit new to this is that a rate of change of T or a final?

Could I get some help with this in terms of getting the script to produce realistic values?

P.S. I have a warning at this time showing up in the output window that says:


2024-05-27 13:21:44.456 (   0.420s) [main            ]           XDMFFile.cpp:250   WARN| No mesh found at '/Xdmf/Domain/Grid[@GridType='Uniform'][1]'. Write mesh before function!

Not quite sure what is meant by it. So far the output.xdmf seems to be crashing the latest desktop version of ParaView if maybe that warning is indicative to some type of problem with the outputting of xdmf. OK… there is also an h5, one thing I might be interested also to know is the meaning of the additional h5 and to its use and contents

2

You haven’t written the mesh to the XDMFFile. This could be dong right after

You are setting your initial condition to Zero

Is this by design of your problem?

Please also look at the solution of your problem visually.

1 Like
3

time t=0:

timet0
timet01601×906 66.4 KB

time t=1:

timet1
timet11328×895 53.8 KB

I had to rescale the color map so the legend will show properly at t=1. It seems to imply so far that things are heating up in places and getting colder in other places.

Yes that is true about the initial condition. It had seemed to me a good enough value at the time to start with just to get the system moving and outputting some solution.

So it seems I have a final temperature and the rate of change of temperature is someplace there in the variational formulation F I would guess so far.

The negative values had thrown me a bit. I feel unsure yet if I should accept such a solution as representative of reality. Does this solution look accurate or is there maybe someplace where I went wrong?

4

Start with a smaller delta t, and only look at the evolution for one time step. You could even reduce this to a 2D problem to begin with.
As you are heating up the boundary with a linear function, there will be a large jump from your initial value (0) to the boundary value (0.1), which might cause some instability.

1 Like
5

OK… One thing I just noticed is that the Dirichlet seems to be set to a function of t and that is not yet what I want. I was thinking ideally the Dirichlet should interpolate between the outside temperature and the temperature dolfinx would like to set at the nodes. I tried my best so far to do that but was all thumbs at it so far…

6

It is unclear what you mean by

Could you write down the mathematical formulation of your bc?

7

This would be the formula:


T_new = 0.5 * (u_prev + T_outside)

So far one of my attempts to compute that would be:


def boundary_condition(u_n, x, t, T_outside=0):
    # Evaluate the previous time step's temperature at the boundary nodes
    u_prev = u_n.eval(x) <---
    # Interpolate between the previous temperature and the outside temperature
    T_new = 0.5 * (u_prev + T_outside)
    return T_new
8

Use a dolfinx.fem.Expression to define this function for interpolation, i.e.

update_expr = dolfinx.fem.Expression(0.5 * (u_prev + T_outside), V.element.interpolation_points())
u_D.interpolate(update_expr)

You can define update_expr outside the temporal loop, and call u_D.interpolate(update_expr) inside the time loop (as it takes into account changes in u_prev (if u_prev is u_n and is a dolfinx.fem.Function).

1 Like