Skip to content
Snippets Groups Projects
Commit d86745d1 authored by Joe Pater's avatar Joe Pater
Browse files

2D case working

parent cfb86c86
No related branches found
No related tags found
No related merge requests found
...@@ -2,10 +2,213 @@ from scipy import integrate as sp_int ...@@ -2,10 +2,213 @@ from scipy import integrate as sp_int
from math import exp from math import exp
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
import numpy as np import numpy as np
from numpy.linalg import inv
import math import math
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
def sum_to(k):
return k*(k+1)/2
class Mesh2D: class Mesh2D:
def __init__(self, R, h, delta_r, delta_h): def __init__(self, R, h, N):
self.delta_r = R / N
self.delta_z = h / N
self.N = int(N)
self.R = R
self.h = h
# Get the index into the global state vector of a particular
# r and z (r = j*delta_r, z = k*delta_z)
def get_idx(self, j, k):
return round((sum_to(self.N) - sum_to(self.N - k)) + j)
def get_state_len(self):
return round(sum_to(self.N))
def is_sloped_boundary(self, j, k):
return (j + k + 1) == self.N
def is_centre_boundary(self, j, k):
return j == 0
def is_bottom_boundary(self, j, k):
return k == 0
def __getitem__(self, x):
return self.get_idx(x[0], x[1])
def get_plot_matrix(self, y):
P = -np.ones((self.N, self.N * 2 - 1))
for k in range(self.N):
for j in range(self.N - k):
P[k, self.N+j-1] = P[k, self.N-j-1] = y[self[j,k]]
return P
def construct_matrix(mesh):
L = mesh.get_state_len()
M = np.zeros((L,L))
dr = mesh.delta_r
dz = mesh.delta_z
for k in range(mesh.N):
for j in range(mesh.N - k):
row = mesh[j,k]
M_up = 1/dz**2
M_down = 1/dz**2
if j != 0:
M_left = 1/dr**2 + 1/(2*j*dr**2)
M_right = 1/dr**2 - 1/(2*j*dr**2)
M_c = -2/dr**2 - 2/dz**2
else:
M_left = 2/dr**2
M_right = 2/dr**2
M_c = -4/dr**2 - 2/dz**2
M[row, mesh[j,k]] = M_c
if j > 0:
M[row, mesh[j-1,k]] = M_right
if k > 0:
M[row, mesh[j,k-1]] = M_down
if not mesh.is_sloped_boundary(j, k):
if mesh.is_centre_boundary(j, k):
M[row, mesh[j+1,k]] = M_left + M_right # Symmetry BC
else:
M[row, mesh[j+1,k]] = M_left
if mesh.is_bottom_boundary(j, k):
M[row, mesh[j,k+1]] = M_up + M_down # Symmetry BC
else:
M[row, mesh[j,k+1]] = M_up
return M
def solve(M, y0, T, delta_t):
#M_inv = inv(delta_t*M + np.eye(len(y0)))
step_mat = inv(np.eye(len(y0)) - delta_t*M)
N_t = round(T / delta_t)
y = y0
sln = np.zeros((N_t, len(y0)))
for i in range(N_t):
if max(y) > 20:
break
y = np.matmul(step_mat, y + delta_t*np.exp(y))
sln[i,:] = y
return sln, i+1
if __name__ == "__main__":
mesh = Mesh2D(50, 50, 40)
# for k in range(mesh.N):
# for j in range(mesh.N-k):
# print("{:<5}".format(mesh[j,k]), end="")
# print("")
Da = 0.001
M = construct_matrix(mesh) / Da
ic_low = 0 #3.4
ic_high = 10 #3.4
i = 0
while i < 100:
i += 1
ic = (ic_low + ic_high) / 2
initial_conditions = ic * np.ones(mesh.get_state_len())
simulation_time = 10
time_step = 0.5
sln, N_t_in_sln = solve(M, initial_conditions, simulation_time, time_step)
print(ic, N_t_in_sln)
if N_t_in_sln < round(simulation_time / time_step):
ic_high = ic
else:
if i > 50:
break
ic_low = ic
print(ic)
num_subplots_target = 8
y_final = sln[-1,:]
P = mesh.get_plot_matrix(y_final)
temp_min = 0.0
temp_max = 6
# Determine subplot grid
ncols = 4
nrows = int(np.ceil(num_subplots_target / ncols)) # For 8 plots, this will be 2x4
fig, axes = plt.subplots(nrows=nrows, ncols=ncols,
figsize=(5 * ncols, 4 * nrows), constrained_layout=True)
axes = axes.ravel() # Flatten for easy indexing
# Choose a colormap (e.g., 'viridis', 'inferno', 'hot', 'coolwarm')
cmap = plt.cm.viridis
# Set the color for values in P that are below temp_min.
# Since your background in P is -1.0, and if temp_min is 0.0,
# these -1.0 values will be rendered as black.
cmap.set_under('white')
# Plot 1: Initial Condition (y0)
ax_current = axes[0]
P_initial = mesh.get_plot_matrix(initial_conditions)
ax_current.set_facecolor('black')
# Store the mappable object from an imshow call for the colorbar
mappable_for_colorbar = ax_current.imshow(P_initial, cmap=cmap, vmin=temp_min,
vmax=temp_max, origin='lower', aspect='auto',
interpolation='nearest')
ax_current.set_title('T = 0.00s (Initial)')
ax_current.set_xticks([])
ax_current.set_yticks([])
plots_created = 1
# Plots 2 to num_subplots_target: Sampled from sln
# We need num_subplots_target - 1 plots from sln
num_plots_from_sln_needed = num_subplots_target - 1
if N_t_in_sln > 0 and num_plots_from_sln_needed > 0:
# Select indices from sln. If N_t_in_sln is less than needed, take all available.
num_to_sample_from_sln = min(num_plots_from_sln_needed, N_t_in_sln)
plot_indices_in_sln = np.linspace(0, N_t_in_sln - 1, num_to_sample_from_sln, dtype=int)
for i in range(num_to_sample_from_sln):
ax_current = axes[plots_created] # Next available subplot
snapshot_idx_in_sln = plot_indices_in_sln[i]
y_snapshot = sln[snapshot_idx_in_sln, :]
P_snapshot = mesh.get_plot_matrix(y_snapshot)
# Time for sln[k] is (k+1)*delta_t
current_time = (snapshot_idx_in_sln + 1) * time_step
ax_current.set_facecolor('black')
ax_current.imshow(P_snapshot, cmap=cmap, vmin=temp_min, vmax=temp_max,
origin='lower', aspect='auto', interpolation='nearest')
ax_current.set_title('T = {:.2f}s'.format(current_time))
ax_current.set_xticks([])
ax_current.set_yticks([])
plots_created += 1
# Hide any remaining unused subplots
for i in range(plots_created, nrows * ncols):
fig.delaxes(axes[i])
# Add a single colorbar for the entire figure
if plots_created > 0:
active_axes = [axes[j] for j in range(plots_created)]
if active_axes:
fig.colorbar(mappable_for_colorbar, ax=active_axes, shrink=0.7, aspect=15*nrows, pad=0.05, label='Temperature')
fig.suptitle('Temperature Distribution Over Time', fontsize=16)
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment