Packaged code by kimmerfx

.gitignore

@@ -0,0 +1,3 @@
+dist/
+kimmerfx_diffusion2d.egg-info/
+__pycache__

README.md

@@ -7,10 +7,29 @@ Please follow the instructions in [pypi_exercise.md](https://github.com/Simulati
 The code used in this exercise is based on [Chapter 7 of the book "Learning Scientific Programming with Python"](https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/).
 
 ## Description
+This code solves the diffusion equation over a two dimensional square domain which is at a certain temperature, and a circular disc at its center which is at a higher temperature. The diffusion equation is solved using the [finite-difference](https://en.wikipedia.org/wiki/Finite_difference_method) method. The thermal diffusivity and initial conditions of the system can be changed by the user. The code produces four plots at various timepoints of the simulation. The diffusion process can be clearly observed in these plots.
 
 ## Installing the package
+```bash
+pip install \
+  --index-url https://test.pypi.org/simple/ \
+  --extra-index-url https://pypi.org/simple \
+  kimmerfx_diffusion2d
+```
 
 ## Running this package
+You can simply run:
+```bash
+diffusion-solve
+```
+in your terminal.
+
+Or use this in a script or in the repl.
+```python
+from kimmerfx_diffusion2d.diffusion2d import solve
+solve()
+```
 
 ## Citing
 
+[Original source (E7.24: The two-dimensional diffusion equation)](https://scipython.com/books/book2/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/)

kimmerfx_diffusion2d/diffusion2d.py

@@ -0,0 +1,84 @@
+"""
+Solving the two-dimensional diffusion equation
+
+Example acquired from https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from .output import create_plot, output_plots
+
+
+def solve(dx: float = 0.1, dy: float = 0.1, D: float = 4.0) -> None:
+    """
+    Solve the 2D diffusion equation and plots the ouput at different simulation times.
+
+    Args:
+        dx (float): Distance between two points in the x direction (mm).
+        dy (float): Distance between two points in the y direction (mm).
+        D (float): Thermal diffusivity of steel, in mm^2/s.
+    """
+    # plate size, mm
+    w = h = 10.0
+
+    # Initial cold temperature of square domain
+    T_cold = 300
+
+    # Initial hot temperature of circular disc at the center
+    T_hot = 1000
+
+    # Number of discrete mesh points in X and Y directions
+    nx, ny = int(w / dx), int(h / dy)
+
+    # Computing a stable time step
+    dx2, dy2 = dx * dx, dy * dy
+    dt = dx2 * dy2 / (2 * D * (dx2 + dy2))
+
+    print("dt = {}".format(dt))
+
+    u0 = T_cold * np.ones((nx, ny))
+    u = u0.copy()
+
+    # Initial conditions - circle of radius r centred at (cx,cy) (mm)
+    r = min(h, w) / 4.0
+    cx = w / 2.0
+    cy = h / 2.0
+    r2 = r**2
+    for i in range(nx):
+        for j in range(ny):
+            p2 = (i * dx - cx) ** 2 + (j * dy - cy) ** 2
+            if p2 < r2:
+                u0[i, j] = T_hot
+
+    def do_timestep(u_nm1, u, D, dt, dx2, dy2):
+        # Propagate with forward-difference in time, central-difference in space
+        u[1:-1, 1:-1] = u_nm1[1:-1, 1:-1] + D * dt * (
+            (u_nm1[2:, 1:-1] - 2 * u_nm1[1:-1, 1:-1] + u_nm1[:-2, 1:-1]) / dx2 + (u_nm1[1:-1, 2:] - 2 * u_nm1[1:-1, 1:-1] + u_nm1[1:-1, :-2]) / dy2
+        )
+
+        u_nm1 = u.copy()
+        return u_nm1, u
+
+    # Number of timesteps
+    nsteps = 101
+    # Output 4 figures at these timesteps
+    n_output = [0, 10, 50, 100]
+    fig_counter = 0
+    fig = plt.figure()
+
+    # Time loop
+    for n in range(nsteps):
+        u0, u = do_timestep(u0, u, D, dt, dx2, dy2)
+
+        # Create figure
+        if n in n_output:
+            fig_counter += 1
+            im = create_plot(fig, fig_counter, T_cold, T_hot, dt, u, n)
+
+    # Plot output figures
+    output_plots(fig, im)
+
+
+if __name__ == "__main__":
+    solve()

kimmerfx_diffusion2d/output.py

@@ -0,0 +1,53 @@
+"""
+Helper functions for diffusion2d.
+
+Create and diplay the plots
+"""
+
+from matplotlib.image import AxesImage
+import matplotlib.pyplot as plt
+from matplotlib.figure import Figure
+import numpy as np
+from numpy.typing import NDArray
+
+
+def create_plot(fig: Figure, fig_counter: int, T_cold: int, T_hot: int, dt: float, u: NDArray[np.floating], n: int) -> AxesImage:
+    """
+    Create a subplot with the current temperature field and title it by time.
+
+    Args:
+        fig (Figure): Matplotlib figure to which the subplot is added.
+        fig_counter (int): Subplot index used with a 2x2 grid layout
+        T_cold (int): Minimum temperature for color scaling (K).
+        T_hot (int): Maximum temperature for color scaling (K).
+        dt (float): Time step size (s).
+        u (NDArray[np.floating]): Temperature field array of shape (nx, ny).
+        n (int): Current time step index.
+
+    Returns:
+        AxesImage: The image object created by `imshow`, suitable for a
+        shared colorbar.
+    """
+    ax = fig.add_subplot(220 + fig_counter)
+    im = ax.imshow(u.copy(), cmap=plt.get_cmap("hot"), vmin=T_cold, vmax=T_hot)  # image for color bar axes
+    ax.set_axis_off()
+    ax.set_title("{:.1f} ms".format(n * dt * 1000))
+    return im
+
+
+def output_plots(fig: Figure, im: AxesImage) -> None:
+    """
+    Finalize the layout, add a colorbar for the last image, and show the figure.
+
+    Args:
+        fig (Figure): Matplotlib figure containing the subplots.
+        im (AxesImage): Image handle used to derive the colorbar scaling.
+
+    Returns:
+        None
+    """
+    fig.subplots_adjust(right=0.85)
+    cbar_ax = fig.add_axes([0.9, 0.15, 0.03, 0.7])
+    cbar_ax.set_xlabel("$T$ / K", labelpad=20)
+    fig.colorbar(im, cax=cbar_ax)
+    plt.show()

pyproject.toml

@@ -0,0 +1,28 @@
+[build-system]
+requires = ["setuptools >= 77.0.3"]
+build-backend = "setuptools.build_meta"
+
+[project]
+version = "0.0.3"
+name="kimmerfx_diffusion2d"
+dependencies = [
+  "numpy",
+  "matplotlib"
+]
+requires-python = ">= 3.8"
+description = "Solves the diffusion equation."
+readme = "README.md"
+license-files = ["LICENSE"]
+keywords = ["diffusion", "2D", "solver"]
+
+classifiers = [
+  "Development Status :: 3 - Alpha",
+]
+
+authors = [
+  {name = "Ishaan Desai"},
+  {name = "Felix Kimmerle"},
+]
+
+[project.scripts]
+diffusion-solve = "kimmerfx_diffusion2d.diffusion2d:solve"