Help with creating own module

[sharanroongta]

I was just trying to create a simple heat conduction example for myself. For 1D, I was able to use the software.

Can you help me with using this for 2D, with an isotropic thermal conductivity matrix? Is it possible? Would I have to change something in the BaseSystem class?

Thanks

[vidarsko]

Sure! By isotropic thermal conductivity, are you referring to the following PDE:

$$ \partial_t T = \alpha \nabla^2 T $$

where $t$ is time, $T$ is temperature and $\alpha$ the thermal diffusivity?

[sharanroongta]

Yes, but in case of 2D, $\alpha$ should be diagonal matrix (tensor) with same value in x-x and y-y direction.

In case of anisotropic case, those 2 values could be different.

[vidarsko]

Let me see if I understand you. Given a thermal conducitvity $\kappa$ (matrix) so that the heat flux $\vec q$ is given by

$$ \vec q = - \kappa \cdot \nabla T. $$

If I take $\kappa$ to be a diagonal matrix matrix with only one component $k$

$$ \kappa = \begin{pmatrix} k & 0 \ 0 & k \ \end{pmatrix}, $$

then the equation reduces to

$$ \vec q = - k \nabla T. $$

The change in temperature $T$ is given by

$$ \partial_t T = - \frac{1}{\rho c_p}\nabla \cdot \vec q, $$

where $\rho$ is the material density and $c_p$ the specific heat, which if taken to be constant gives

$$ \partial_t T = \frac{k}{\rho c_p} \nabla^2 T \equiv \alpha \nabla^2 T, $$

meaning that if you have isotropic heat conduction in a uniform media, the heat equation reduces to the above, where $\alpha$ is a number, not a matrix no?

Or do you want a computational setup where you can modify the matrix $\kappa$ to be possibly non-isotropic? Both are possible.

[sharanroongta]

Or do you want a computational setup where you can modify the matrix 
 to be possibly non-isotropic? Both are possible.

Exactly what you said in the last sentence. The flexibility to simulate anisotropic conditions. Can you guide me how I can do that?

[vidarsko]

That is a very nice problem, which we need to formulate in terms of a PDE. The thermal conductivity is a symmetric matrix with three independent components

$$ \kappa = \begin{pmatrix} \kappa{11} & \kappa{12} \ \kappa{12} & \kappa{22}\ \end{pmatrix} $$

The heat flux becomes

$$ \begin{pmatrix} q_1 \ q2 \end{pmatrix} =\begin{pmatrix} \kappa{11} & \kappa{12} \ \kappa{12} & \kappa_{22}\ \end{pmatrix} \begin{pmatrix} \partial_x T \ \partialy T \end{pmatrix} =\begin{pmatrix} \kappa{11} \partialx T + \kappa{12} \partialy T \ \kappa{12} \partialx T + \kappa{22} \partial_y T \ \end{pmatrix} $$

Assuming a uniform media with uniform heat capacity, we get

$$ \partial_t T = \frac{1}{\rho c_p} \nabla \cdot \vec q = \frac{1}{\rho cp} (\kappa{11} \partial{xx} T + 2 \kappa{12} \partial{xy} T + \kappa{22} \partial_{yy} T ) $$

If we define the thermal diffusivity matrix as $\alpha_{ij} = \frac{1}{\rho cp} \kappa{ij}$, we get

$$ \partialt T = \alpha{11} \partial{xx} T + 2 \alpha{12} \partial{xy} T + \alpha{22} \partial_{yy} T $$

To create a system with comfit implementing this equation, you create a class like in the "make your own model" tutorial.

class HeatEquationSystem(cf.BaseSystem):

  def __init__(self, dim, alpha11, alpha12, alpha22, **kwargs):
    self.alpha11 = alpha11
    self.alpha12 = alpha12
    self.alpha22 = alpha22
    super().__init__(dim, **kwargs)

This equation is purely linear, so you define the linear and non linear parts

def calc_omega_f(self):
    return self.alpha11*self.dif[0]**2 + 2*self.alpha12*self.dif[0]*self.dif[1] + self.alpha22*self.dif[1]**2

# Add method to class
HeatEquationSystem.calc_omega_f = calc_omega_f

# Make the non-linear operator
def calc_nonlinear_evolution_function_f(self, field, t):
    return 0

# Add method to class
HeatEquationSystem.calc_nonlinear_evolution_function_f = calc_nonlinear_evolution_function_f

Note: I implemented this to be specific to 2 dimensions, specifying alpha11, alpha12 etc. It is possible to write it more generally using a summation over components which would make it more easily extendible to 3 dimensions. Now, the solver

import numpy as np

def evolve(self, number_steps):
    omega_f = calc_omega_f(self)

    integrating_factors_f, solver = self.calc_integrating_factors_f_and_solver(omega_f, 'ETD2RK')

    for n in range(number_steps):
        self.T, self.T_f = solver(integrating_factors_f, 
                                    self.calc_nonlinear_evolution_function_f, 
                                    self.T, self.T_f)
        self.T = np.real(self.T)

# Add evolve method to class
HeatEquationSystem.evolve = evolve

Now you can instanciate a system and set an initial condition

hes = HeatEquationSystem(2,1,0,0.01)

import scipy as sp

hes.T = hes.calc_Gaussian()
hes.T_f = sp.fft.fftn(hes.T)

Plot the intial condition

hes.plot_field(hes.T)

Evolve

hes.evolve(1000)

And plot, to see htat in this case, there is almost not change in the temperature in the $y$-direction

hes.plot_field(hes.T)

I was able to run this code in colab. Let me know if you have further questions!

[sharanroongta]

Super! Now I understand, many thanks. I would just give it a run once.