Solver for heat equation with a source term and different boundary conditions

[abhishek-appana]

Hello everyone,

I am looking to solve the heat equation on a rectangular geometry with different boundary conditions on the edges as explained below:

The heat equation is:
$\rho C_p \partial T/\partial t = k \Delta T + Q_{gen} $

The edges 1, 2, 3, 5, 7, 8 have the following boundary condition:
$-k dT/dn = h(T-T_{inf})$
Edges 4 & 6 have the boundary condition:
$-k dT/dn = q_{tabs}$

Below is the code I have written for this. The figure gets generated for t=0, but the solution doesn't get published/generated from the next time step. Can someone please let me know what needs to be changed?

Thank you

from typing import Iterator, Tuple

import numpy as np
from scipy.sparse.linalg import splu

from skfem import *
from skfem.helpers import dot, grad

length = 0.5
height = 0.5
ncells = 50

T_inf = 298
Q_gen = 1000
q_tabs = 500
rho = 25
C_p = 90
k = 25
h = 10

dt = 0.01

mesh = (MeshQuad
        .init_tensor(np.linspace(0, length, ncells),
                     np.linspace(0, height, ncells))
        .with_boundaries({'1': lambda x: np.isclose(x[1],0),
                          '2': lambda x: np.isclose(x[0],0),
                          '3': lambda x: np.isclose(x[1],0.5) & np.greater_equal(x[0],0) & np.less_equal(x[0],0.1),
                          '4': lambda x: np.isclose(x[1],0.5) & np.greater_equal(x[0],0.1) & np.less_equal(x[0],0.2),
                          '5': lambda x: np.isclose(x[1],0.5) & np.greater_equal(x[0],0.2) & np.less_equal(x[0],0.3),
                          '6': lambda x: np.isclose(x[1],0.5) & np.greater_equal(x[0],0.3) & np.less_equal(x[0],0.4),
                          '7': lambda x: np.isclose(x[1],0.5) & np.greater_equal(x[0],0.4) & np.less_equal(x[0],0.5),
                          '8': lambda x: np.isclose(x[0],0.5)}))

e = ElementQuad2()  # or ElementQuad1
basis = {
    'cell': Basis(mesh, e),
    'pos_tab': FacetBasis(mesh, e, facets=mesh.boundaries['4']),
    'neg_tab': FacetBasis(mesh, e, facets=mesh.boundaries['6']),
    'bottom': FacetBasis(mesh, e, facets=mesh.boundaries['1']),
    'left': FacetBasis(mesh, e, facets=mesh.boundaries['2']),
    'top_1': FacetBasis(mesh, e, facets=mesh.boundaries['3']),
    'top_2': FacetBasis(mesh, e, facets=mesh.boundaries['5']),
    'top_3': FacetBasis(mesh, e, facets=mesh.boundaries['7']),
    'right': FacetBasis(mesh, e, facets=mesh.boundaries['8'])}

@BilinearForm
def laplace(u, v, _):
    return dot(grad(u),grad(v))

@BilinearForm
def mass(u, v, _):
    return u*v

@LinearForm
def tabs(v, _):
    return dt*q_tabs*v

@BilinearForm
def rest(u, v, _):
    return dt*h*(u - T_inf)*v

@LinearForm
def loss(v, _):
    return dt*Q_gen*v

T_init = np.repeat(298, basis['cell'].doflocs.shape[1])

M = rho * C_p * asm(mass,basis['cell'])
L = k * asm(laplace, basis['cell'])
rhs_1 = -1 * (asm(tabs,basis['neg_tab']) + asm(tabs, basis['pos_tab']))
rhs_2 = (-1 * (asm(rest,basis['bottom']) + asm(rest,basis['left']) + asm(rest,basis['top_1']) 
              + asm(rest,basis['top_2']) + asm(rest,basis['top_3']) + asm(rest,basis['right'])))
rhs_3 = asm(loss, basis['cell'])

A = M + 0.5 * L * dt
B = M - 0.5 * L * dt

backsolve = splu(A.T).solve

def evolve(t: float, 
           u: np.ndarray) -> Iterator[Tuple[float, np.ndarray]]:
    
    while t < 10:
        t, u = t + dt, backsolve((B @ u) + rhs_1 + rhs_2 + rhs_3)
        yield t, u


if __name__ == '__main__':

    from argparse import ArgumentParser
    from pathlib import Path

    from matplotlib.animation import FuncAnimation
    import matplotlib.pyplot as plt

    from skfem.visuals.matplotlib import plot

    parser = ArgumentParser(description='heat equation in a rectangle')
    parser.add_argument('-g', '--gif', action='store_true', 
                        help='write animated GIF', )
    args = parser.parse_args()

    ax = plot(mesh, T_init[basis['cell'].nodal_dofs.flatten()], shading='gouraud')
    title = ax.set_title('t = 0.00')
    field = ax.get_children()[0]  # vertex-based temperature-colour
    fig = ax.get_figure()
    fig.colorbar(field)


    def update(event):
        t, u = event

        title.set_text(f'$t$ = {t:.2f}')
        field.set_array(u[basis['cell'].nodal_dofs.flatten()])

    animation = FuncAnimation(
        fig,
        update,
        evolve(0., T_init),
        repeat=False,
        interval=50,
    )
    if args.gif:
        animation.save(Path(__file__).with_suffix('.gif'), 'imagemagick')
    else:
        plt.show()

[kinnala]

This part does not make sense:

@BilinearForm
def rest(u, v, _):
    return dt*h*(u - T_inf)*v

It is clearly not bilinear, i.e., linear both in $u$ and $v$.
I suggest you carefully write down the weak formulation of your problem.

[gdmcbain]

Yes, switching from Crank-Nicolson (theta=0.5) to backward Euler (theta=1) is a good idea; much more stable, much less sensitive to initial conditions. No need to rewrite the code (from e.g. ex19) just change the value of theta.

[kinnala]

I'd probably try reproducing the static case $$\partial_t T = 0$$ which should correspond to $$t \rightarrow \infty$$ because you have no time-dependent loading or BC, using the time-dependent solver by running it long enough, just to see if the time discretization is correct.

[gdmcbain]

I'd probably try reproducing the static case $$\partial_t T = 0$$ which should correspond to $$t \rightarrow \infty$$ because you have no time-dependent loading or BC, using the time-dependent solver by running it long enough, just to see if the time discretization is correct.

Yes, with theta=1 ("implicit" or "backward" Euler), the time step dt may been taken as long as one likes, even to numpy.Inf if one divides through by dt, in which case the one backward Euler step reduces to the steady state solve.

[gdmcbain]

The actual problem I am trying to solve is slightly more complicated than this. (…term is not uniform across the geometry and the geometry is 3 dimensional). I am not sure what else I can do to simplify the example.

Another idea here is to avoid defining the .mesh_boundaries via numpy.isclose. This tends to be complicated and a source of errors. A tidier way is to define the geometry and meshes in a specialized tool like Gmsh, including Physical Curve entities for the different parts of the mixed boundary condition, export in Gmsh's native MSH file format, and load with Mesh.load as in ex28:

scikit-fem/docs/examples/ex28.py

Line 72 in 2bf93b7

mesh = Mesh.load(Path(__file__).parent / 'meshes' / 'ex28.msh')

Defining the geometry and mesh in something like Gmsh allows the geometry to be defined exactly, including the patches, and then meshed conformally; this avoids having to try to define the boundary patches quasi-geometrically on a discrete mesh without being able to refer to where the nodes and elements of the mesh actually area.

This will become even more important in three dimensions.

[abhishek-appana]

Thank you for your suggestions @kinnala @gdmcbain.

I am exploring all these and some other ideas I have. Will update the forum!

[kinnala]

@kinnala: Thank you for your suggestion. The results from this are not at all close to what I expected. I will look into using this method further and see how it goes!

I have never carefully studied time-dependent PDE problems other than learning the basics so I cannot help too much in the implementation or verification. Above I've attempted to use the implicit Euler method which I learned as a student long time ago but haven't used since. It is very simple but I did not verify the result in any way. There must be some test problems available for the heat equation where you can verify that you have implemented the solver correctly, unfortunately I have no time to dig through the references right now.

I must say that the package scikit-fem is primarily for the spatial discretization, which I am more familiar with, even though others have demonstrated the application of scikit-fem to the spatial discretization of a time-dependent PDE problem. I'm quite sure many other integrated solvers exist for time-dependent problems that have thoroughly tested physics modules for the heat equation.

[gdmcbain]

Yes examples 19 & 39 demonstrate the generalized trapezoidal rule for the heat equation in one and two dimensions against the exact solutions, as derived using the classical techniques of Carslaw & Jaeger's Conduction of Heat in Solids; no problem there.