Implementation of heat conduction equation with robin and Neumann boundary conditions

[abhishek-appana]

Hello everyone,

I am working to solve the heat conduction equation with Robin and Neumann boundary conditions on different edges of the boundary. I am trying to modify the code given for example 19 to make it happen. I am not able to understand some parts of that code (lines 62 & 67).

What does 'mass' represent on line 62?

What is the function 'penalize' on line 67?

Can someone direct me to any available documentation on these?

Thank you!

[gdmcbain]

mass is the "mass" linear form which when assembled for a given basis gives the Gramian of that basis; here M. This M is usually called the "mass matrix" in the finite element method; e.g. Ern & Guermond (2004, Theory and Practice of Finite Elements; §3.3.4). It invariably arises in transient or eigenvalue problems. Definition:

scikit-fem/skfem/models/poisson.py

Line 18 in 68bd6fe

def mass(u, v, _):

penalize imposes the essential boundary condition, here zero, by dividing excursions by a very small number to produce a very big residual unless the value on the boundary is almost the desired value. It's an alternative to enforce and condense. Penalization is discussed by E. & G. in §8.4.4.

[gdmcbain]

The last term looks like a volumetric source. This is handled as in

scikit-fem/docs/examples/ex01.py

Lines 23 to 26 in 6159aa3

	# this method could also be imported from skfem.models.unit_load
	@LinearForm
	def rhs(v, _):
	return 1.0 * v

(the source being unit in ex01).

The penultimate term is a superficial source; an example can be found in:

scikit-fem/docs/examples/ex28.py

Line 102 in 6159aa3

b = (asm(unit_load, basis['heated'])

(again for a unit load, but the LinearForm can be modified for nonuniform sources).

[abhishek-appana]

Hello @gdmcbain,

Thank you for your response! That wasn't my question.

I understand how to express them as linear/bilinear forms and assemble the terms. My question was in terms of solving the problem.

Below is the code I have written for this. It is not working. I am not sure what is wrong!

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

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

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 = asm(loss, basis['cell'])

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

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)
        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()

[gdmcbain]

Right, yes, sorry. But when you say "it is not working", how do you mean? What happens? I can't run this here because a number of terms are undefined (rho, C_p, basis, …).

[abhishek-appana]

The result is displayed for t=0. But it is not getting solved post that.

Sorry about that. Here is the full code. It's slightly different from above (Had to add another rhs term)

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()

[abhishek-appana]

I will start a new discussion about this! Helps track it better