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

Incorporated 1D case into 2D code

parent cd036439
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
diagrams.drawio
+ 103
− 1
View file @ 42ec5bc6
<mxfile host="Electron" agent="Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/26.0.16 Chrome/132.0.6834.196 Electron/34.2.0 Safari/537.36" version="26.0.16"> <mxfile host="Electron" agent="Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) draw.io/26.0.16 Chrome/132.0.6834.196 Electron/34.2.0 Safari/537.36" version="26.0.16">
<diagram name="Page-1" id="ns3onrHNIQr2ljXGaBcB"> <diagram name="Page-1" id="ns3onrHNIQr2ljXGaBcB">
<mxGraphModel dx="367" dy="217" grid="0" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="3300" pageHeight="4681" math="0" shadow="0"> <mxGraphModel dx="392" dy="232" grid="0" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="3300" pageHeight="4681" math="0" shadow="0">
<root> <root>
<mxCell id="0" /> <mxCell id="0" />
<mxCell id="1" parent="0" /> <mxCell id="1" parent="0" />
...@@ -480,6 +480,108 @@ ...@@ -480,6 +480,108 @@
<mxCell id="1nF5tneqOkAijIU8bpes-103" value="s&lt;sub&gt;i+1&lt;/sub&gt;" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;fontSize=8;" parent="1" vertex="1"> <mxCell id="1nF5tneqOkAijIU8bpes-103" value="s&lt;sub&gt;i+1&lt;/sub&gt;" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;fontSize=8;" parent="1" vertex="1">
<mxGeometry x="936" y="804" width="9" height="10" as="geometry" /> <mxGeometry x="936" y="804" width="9" height="10" as="geometry" />
</mxCell> </mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-1" value="5" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FF66FF;" vertex="1" parent="1">
<mxGeometry x="490" y="970" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-2" value="4" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FF9999;" vertex="1" parent="1">
<mxGeometry x="470" y="970" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-3" value="3" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FF9999;" vertex="1" parent="1">
<mxGeometry x="450" y="970" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-4" value="2" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FF9999;" vertex="1" parent="1">
<mxGeometry x="430" y="970" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-5" value="1" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#66FFFF;" vertex="1" parent="1">
<mxGeometry x="410" y="970" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-6" value="9" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#6666FF;" vertex="1" parent="1">
<mxGeometry x="470" y="950" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-7" value="8" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#66FF66;" vertex="1" parent="1">
<mxGeometry x="450" y="950" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-8" value="7" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#66FF66;" vertex="1" parent="1">
<mxGeometry x="430" y="950" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-9" value="6" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FFFF99;" vertex="1" parent="1">
<mxGeometry x="410" y="950" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-10" value="12" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#6666FF;" vertex="1" parent="1">
<mxGeometry x="450" y="930" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-11" value="11" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#66FF66;" vertex="1" parent="1">
<mxGeometry x="430" y="930" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-12" value="10" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FFFF99;" vertex="1" parent="1">
<mxGeometry x="410" y="930" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-13" value="14" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#6666FF;" vertex="1" parent="1">
<mxGeometry x="430" y="910" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-14" value="13" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FFFF99;" vertex="1" parent="1">
<mxGeometry x="410" y="910" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-15" value="15" style="whiteSpace=wrap;html=1;aspect=fixed;fillColor=#FFB366;" vertex="1" parent="1">
<mxGeometry x="410" y="890" width="20" height="20" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-16" value="j = 0" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="380" y="990" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-17" value="j = 1" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="400" y="990" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-18" value="j = 2" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="420" y="990" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-19" value="j = 3" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="440" y="990" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-20" value="j = 4" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="460" y="990" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-21" value="k = 0" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="364" y="970" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-23" value="k = 1" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="364" y="950" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-24" value="k = 2" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="364" y="930" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-25" value="k = 3" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="364" y="910" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-26" value="k = 4" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;rotation=-45;" vertex="1" parent="1">
<mxGeometry x="364" y="890" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-27" value="N = 5" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;fontStyle=5" vertex="1" parent="1">
<mxGeometry x="457" y="890" width="60" height="30" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-29" value="Case 6" style="whiteSpace=wrap;html=1;fillColor=#FFB366;" vertex="1" parent="1">
<mxGeometry x="530" y="987.5" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-31" value="Case 1" style="whiteSpace=wrap;html=1;fillColor=#66FF66;" vertex="1" parent="1">
<mxGeometry x="530" y="900" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-32" value="Case 2" style="whiteSpace=wrap;html=1;fillColor=#6666FF;" vertex="1" parent="1">
<mxGeometry x="530" y="917.5" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-33" value="Case 3" style="whiteSpace=wrap;html=1;fillColor=#FFFF99;" vertex="1" parent="1">
<mxGeometry x="530" y="935" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-34" value="Case 4" style="whiteSpace=wrap;html=1;fillColor=#FF9999;" vertex="1" parent="1">
<mxGeometry x="530" y="952.5" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-35" value="Case 5" style="whiteSpace=wrap;html=1;fillColor=#66FFFF;" vertex="1" parent="1">
<mxGeometry x="530" y="970" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-36" value="Case 7" style="whiteSpace=wrap;html=1;fillColor=#FF66FF;" vertex="1" parent="1">
<mxGeometry x="530" y="1005" width="80" height="15" as="geometry" />
</mxCell>
<mxCell id="l6ZRUSDHVu6_06W502uG-37" value="Key:" style="text;html=1;align=center;verticalAlign=middle;whiteSpace=wrap;rounded=0;fontStyle=1" vertex="1" parent="1">
<mxGeometry x="515" y="875" width="60" height="30" as="geometry" />
</mxCell>
</root> </root>
</mxGraphModel> </mxGraphModel>
</diagram> </diagram>
......
fire2d.py
+ 161
− 21
View file @ 42ec5bc6
...@@ -11,6 +11,21 @@ from scipy.optimize import minimize ...@@ -11,6 +11,21 @@ from scipy.optimize import minimize
def sum_to(k): def sum_to(k):
return k*(k+1)/2 return k*(k+1)/2
class Mesh1D:
def __init__(self, L, N):
self.delta_x = L / N
self.N = int(N)
self.L = L
def get_idx(self, j):
return j
def get_state_len(self):
return self.N
def __getitem__(self, x):
return self.get_idx(x)
class Mesh2D: class Mesh2D:
def __init__(self, R, h, N): def __init__(self, R, h, N):
self.delta_r = R / N self.delta_r = R / N
...@@ -48,7 +63,7 @@ class Mesh2D: ...@@ -48,7 +63,7 @@ class Mesh2D:
return P return P
def construct_matrix(mesh): def discrete_laplacian_2d_cyl(mesh):
L = mesh.get_state_len() L = mesh.get_state_len()
M = np.zeros((L,L)) M = np.zeros((L,L))
dr = mesh.delta_r dr = mesh.delta_r
...@@ -88,6 +103,27 @@ def construct_matrix(mesh): ...@@ -88,6 +103,27 @@ def construct_matrix(mesh):
return M return M
def discrete_laplacian_1d(mesh):
L = mesh.get_state_len()
M = np.zeros((L,L))
dx = mesh.delta_x
for j in range(mesh.N):
row = mesh[j]
M_left = 1/dx**2
M_c = -2/dx**2
M_right = 1/dx**2
M[row, mesh[j]] = M_c
if j != 0:
M[row, mesh[j-1]] = M_left
if j != mesh.N-1:
M[row, mesh[j+1]] = M_right
return M
def solve(M, y0, T, delta_t): def solve(M, y0, T, delta_t):
#M_inv = inv(delta_t*M + np.eye(len(y0))) #M_inv = inv(delta_t*M + np.eye(len(y0)))
step_mat = inv(np.eye(len(y0)) - delta_t*M) step_mat = inv(np.eye(len(y0)) - delta_t*M)
...@@ -118,38 +154,40 @@ def solve(M, y0, T, delta_t): ...@@ -118,38 +154,40 @@ def solve(M, y0, T, delta_t):
# print(ymin) # print(ymin)
# return ymin.x # return ymin.x
if __name__ == "__main__": def run_2d():
mesh = Mesh2D(50, 50, 60) mesh = Mesh2D(1, 1, 50)
# for k in range(mesh.N): # for k in range(mesh.N):
# for j in range(mesh.N-k): # for j in range(mesh.N-k):
# print("{:<5}".format(mesh[j,k]), end="") # print("{:<5}".format(mesh[j,k]), end="")
# print("") # print("")
Da = 0.001 ic = 0
M = construct_matrix(mesh) / Da
ic_low = 0 #3.4
ic_high = 0 #3.4 Da_low = 0
Da_high = 10
i = 0 i = 0
while i < 100: while i < 100:
i += 1 i += 1
ic = (ic_low + ic_high) / 2 Da = (Da_low + Da_high) / 2
M = discrete_laplacian_2d_cyl(mesh) / Da
initial_conditions = ic * np.ones(mesh.get_state_len()) initial_conditions = ic * np.ones(mesh.get_state_len())
simulation_time = 10 simulation_time = 100
time_step = 0.5 time_step = 0.01
sln, N_t_in_sln = solve(M, initial_conditions, simulation_time, time_step) sln, N_t_in_sln = solve(M, initial_conditions, simulation_time, time_step)
print(ic, N_t_in_sln) print(Da, N_t_in_sln)
#break
if N_t_in_sln < round(simulation_time / time_step): if N_t_in_sln < round(simulation_time / time_step):
ic_high = ic Da_high = Da
else: else:
if i > 0: if i > 20:
break break
ic_low = ic Da_low = Da
print(ic) print(Da)
num_subplots_target = 8 num_subplots_target = 8
...@@ -157,7 +195,7 @@ if __name__ == "__main__": ...@@ -157,7 +195,7 @@ if __name__ == "__main__":
P = mesh.get_plot_matrix(y_final) P = mesh.get_plot_matrix(y_final)
temp_min = 0.0 temp_min = 0.0
temp_max = 0.5 temp_max = 5
# Determine subplot grid # Determine subplot grid
ncols = 4 ncols = 4
...@@ -183,7 +221,7 @@ if __name__ == "__main__": ...@@ -183,7 +221,7 @@ if __name__ == "__main__":
mappable_for_colorbar = ax_current.imshow(P_initial, cmap=cmap, vmin=temp_min, mappable_for_colorbar = ax_current.imshow(P_initial, cmap=cmap, vmin=temp_min,
vmax=temp_max, origin='lower', aspect='auto', vmax=temp_max, origin='lower', aspect='auto',
interpolation='nearest') interpolation='nearest')
ax_current.set_title('T = 0.00s (Initial)') ax_current.set_title('T = 0.00 (Initial)')
ax_current.set_xticks([]) ax_current.set_xticks([])
ax_current.set_yticks([]) ax_current.set_yticks([])
...@@ -197,7 +235,7 @@ if __name__ == "__main__": ...@@ -197,7 +235,7 @@ if __name__ == "__main__":
# Select indices from sln. If N_t_in_sln is less than needed, take all available. # 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) 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) plot_indices_in_sln = np.linspace(0, N_t_in_sln - 2, num_to_sample_from_sln, dtype=int)
for i in range(num_to_sample_from_sln): for i in range(num_to_sample_from_sln):
ax_current = axes[plots_created] # Next available subplot ax_current = axes[plots_created] # Next available subplot
...@@ -212,7 +250,7 @@ if __name__ == "__main__": ...@@ -212,7 +250,7 @@ if __name__ == "__main__":
ax_current.set_facecolor('black') ax_current.set_facecolor('black')
ax_current.imshow(P_snapshot, cmap=cmap, vmin=temp_min, vmax=temp_max, ax_current.imshow(P_snapshot, cmap=cmap, vmin=temp_min, vmax=temp_max,
origin='lower', aspect='auto', interpolation='nearest') origin='lower', aspect='auto', interpolation='nearest')
ax_current.set_title('T = {:.2f}s'.format(current_time)) ax_current.set_title('T = {:.2f}'.format(current_time))
ax_current.set_xticks([]) ax_current.set_xticks([])
ax_current.set_yticks([]) ax_current.set_yticks([])
plots_created += 1 plots_created += 1
...@@ -227,5 +265,107 @@ if __name__ == "__main__": ...@@ -227,5 +265,107 @@ if __name__ == "__main__":
if active_axes: if active_axes:
fig.colorbar(mappable_for_colorbar, ax=active_axes, shrink=0.7, aspect=15*nrows, pad=0.05, label='Temperature') 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) fig.suptitle('Temperature Distribution Over Time')
plt.show() plt.savefig("1d-temp.pdf", bbox_inches="tight")
print("Done")
#plt.show()
def run_1d():
mesh = Mesh1D(2, 500)
ic = 0
Da_low = 0 #3.4
Da_high = 2 #3.4
i = 0
while i < 100:
i += 1
Da = (Da_low + Da_high) / 2
M = discrete_laplacian_1d(mesh) / Da
initial_conditions = ic * np.ones(mesh.get_state_len())
simulation_time = 100
time_step = 0.01
sln, N_t_in_sln = solve(M, initial_conditions, simulation_time, time_step)
print(Da, N_t_in_sln)
#break
if N_t_in_sln < round(simulation_time / time_step):
Da_high = Da
else:
if i > 20:
break
Da_low = Da
print(ic)
num_subplots_target = 8
y_final = sln[-1,:]
temp_min = 0.0
temp_max = 3
# 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
# Plot 1: Initial Condition (y0)
ax_current = axes[0]
P_initial = initial_conditions
ax_current.plot(initial_conditions)
ax_current.set_title('T = 0.00 (Initial)')
ax_current.set_xticks([])
#ax_current.set_yticks([])
ax_current.set_ylim([temp_min, temp_max])
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, :]
# Time for sln[k] is (k+1)*delta_t
current_time = (snapshot_idx_in_sln + 1) * time_step
ax_current.plot(y_snapshot)
ax_current.set_title('T = {:.2f}'.format(current_time))
ax_current.set_xticks([])
if (plots_created % 4) != 0:
ax_current.set_yticks([])
ax_current.set_ylim([temp_min, temp_max])
plots_created += 1
# Hide any remaining unused subplots
for i in range(plots_created, nrows * ncols):
fig.delaxes(axes[i])
fig.suptitle('Temperature Distribution Over Time')
plt.savefig("1d-temp.pdf", bbox_inches="tight")
#plt.show()
print("Done")
if __name__ == "__main__":
plt.rcParams.update({'font.size': 22})
run_2d()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment