diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..06d6d439 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +__pycache__/ +build/ +dist/ +*.egg-info/ \ No newline at end of file diff --git a/README.md b/README.md index fe548354..d0580f29 100644 --- a/README.md +++ b/README.md @@ -8,12 +8,33 @@ The code used in this exercise is based on [Chapter 7 of the book "Learning Scie ## Project description +This code simulates the 2D diffusion equation on a square domain, with the entire domain set to an initial temperature and a circular disc at the center maintained at a higher temperature. +Using the Finite Difference Method, the code solves the diffusion equation, allowing the user to adjust the thermal diffusivity and initial conditions of the system. +The simulation outputs four plots at different time points, vividly illustrating the progression of the diffusion process. + ## Installing the package ### Using pip3 to install from PyPI +```bash +pip install --user --index-url https://test.pypi.org/simple/ tischlre_diffusion2d +``` + ### Required dependencies +```bash +pip install numpy matplotlib +``` + ## Running this package +```python +from tischlre_diffusion2d.diffusion2d import solve + +solve(dx = 0.1, dy = 0.1, D = 4) +``` + ## Citing + + +This was forked from https://github.com/Simulation-Software-Engineering/diffusion2D. \ No newline at end of file diff --git a/diffusion2d.py b/diffusion2d.py deleted file mode 100644 index c0c6083a..00000000 --- a/diffusion2d.py +++ /dev/null @@ -1,81 +0,0 @@ -""" -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 - -# plate size, mm -w = h = 10. -# intervals in x-, y- directions, mm -dx = dy = 0.1 -# Thermal diffusivity of steel, mm^2/s -D = 4. - -# Initial cold temperature of square domain -T_cold = 300 - -# Initial hot temperature of circular disc at the center -T_hot = 700 - -# 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 - 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)) - -# Plot output figures -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() diff --git a/figures/initial.png b/figures/initial.png new file mode 100644 index 00000000..42b9491d Binary files /dev/null and b/figures/initial.png differ diff --git a/figures/testpypI.png b/figures/testpypI.png new file mode 100644 index 00000000..60c8e1d4 Binary files /dev/null and b/figures/testpypI.png differ diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..9429509f --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,26 @@ +[build-system] +requires = ["setuptools", "wheel"] + +[project] +name = "tischlre_diffusion2d" +version = "0.0.1" +description = "An example implementation for solving diffusion equation in 2D" +readme = "README.md" +requires-python = ">=3.6" +license = { file = "LICENSE" } +keywords = ["sse", "simulation", "diffusion"] +classifiers = [ + "Programming Language :: Python :: 3", + "Operating System :: OS Independent", +] +dependencies = [ + "numpy", + "matplotlib", +] + +[project.urls] +Homepage = "https://github.com/Simulation-Software-Engineering/diffusion2D" +Repository = "https://github.com/Simulation-Software-Engineering/diffusion2D" + +[project.entry-points."simulation"] +solve = "tischlre_diffusion2d:solve" diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..c40e2f99 --- /dev/null +++ b/setup.py @@ -0,0 +1,11 @@ +from setuptools import setup, find_packages + +with open("README.md", "r", encoding="utf-8") as fh: + long_description = fh.read() + +if __name__ == "__main__": + setup( + long_description=long_description, + long_description_content_type="text/markdown", + packages=find_packages(exclude=["figures"]), + ) \ No newline at end of file diff --git a/test.py b/test.py new file mode 100644 index 00000000..b1234474 --- /dev/null +++ b/test.py @@ -0,0 +1,3 @@ +from tischlre_diffusion2d.diffusion2d import solve + +solve(dx = 0.1, dy = 0.1, D = 4) \ No newline at end of file diff --git a/tischlre_diffusion2d/__init__.py b/tischlre_diffusion2d/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tischlre_diffusion2d/diffusion2d.py b/tischlre_diffusion2d/diffusion2d.py new file mode 100644 index 00000000..b0a8e46e --- /dev/null +++ b/tischlre_diffusion2d/diffusion2d.py @@ -0,0 +1,75 @@ +""" +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 tischlre_diffusion2d.output import create_plot, output_plots + +def solve(dx = 0.1, dy = 0.1, D = 4): + # plate size, mm + w = h = 10. + # intervals in x-, y- directions, mm + dx = dy = 0.1 + # Thermal diffusivity of steel, mm^2/s + D = 4. + + # Initial cold temperature of square domain + T_cold = 300 + + # Initial hot temperature of circular disc at the center + T_hot = 700 + + # 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: + im, fig_counter = create_plot(u, T_cold, T_hot, n, dt, fig_counter, fig) + + # Plot output figures + output_plots(fig, im) \ No newline at end of file diff --git a/tischlre_diffusion2d/output.py b/tischlre_diffusion2d/output.py new file mode 100644 index 00000000..46483c92 --- /dev/null +++ b/tischlre_diffusion2d/output.py @@ -0,0 +1,16 @@ +import matplotlib.pyplot as plt + +def create_plot(u, T_cold, T_hot, n, dt, fig_counter, fig): + fig_counter += 1 + ax = fig.add_subplot(220 + fig_counter) + im = ax.imshow(u.copy(), cmap=plt.get_cmap('hot'), vmin=T_cold, vmax=T_hot) + ax.set_axis_off() + ax.set_title('{:.1f} ms'.format(n * dt * 1000)) + return im, fig_counter + +def output_plots(fig, im): + 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() \ No newline at end of file