diff --git a/diagrams.drawio b/diagrams.drawio
index ffc0c9fd9bab570b782631e4478f07febdbd3650..7acfc633f0c907d521958cb9b5002086ad9d2e61 100644
--- a/diagrams.drawio
+++ b/diagrams.drawio
@@ -1,6 +1,6 @@
<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">
- <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>
<mxCell id="0" />
<mxCell id="1" parent="0" />
@@ -480,6 +480,108 @@
<mxCell id="1nF5tneqOkAijIU8bpes-103" value="s<sub>i+1</sub>" 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" />
</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>
</mxGraphModel>
</diagram>
diff --git a/fire2d.py b/fire2d.py
index dafd15d6e24a2798dc9c8ac76b6098b8ff397f58..0364c88eedbeeb3386eb372996c211e1ef7e9aa7 100644
--- a/fire2d.py
+++ b/fire2d.py
@@ -11,6 +11,21 @@ from scipy.optimize import minimize
def sum_to(k):
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:
def __init__(self, R, h, N):
self.delta_r = R / N
@@ -48,7 +63,7 @@ class Mesh2D:
return P
-def construct_matrix(mesh):
+def discrete_laplacian_2d_cyl(mesh):
L = mesh.get_state_len()
M = np.zeros((L,L))
dr = mesh.delta_r
@@ -88,6 +103,27 @@ def construct_matrix(mesh):
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):
#M_inv = inv(delta_t*M + np.eye(len(y0)))
step_mat = inv(np.eye(len(y0)) - delta_t*M)
@@ -118,38 +154,40 @@ def solve(M, y0, T, delta_t):
# print(ymin)
# return ymin.x
-if __name__ == "__main__":
- mesh = Mesh2D(50, 50, 60)
+def run_2d():
+ mesh = Mesh2D(1, 1, 50)
# 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 = 0
+
- ic_low = 0 #3.4
- ic_high = 0 #3.4
+ Da_low = 0
+ Da_high = 10
i = 0
while i < 100:
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())
- simulation_time = 10
- time_step = 0.5
+ simulation_time = 100
+ time_step = 0.01
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):
- ic_high = ic
+ Da_high = Da
else:
- if i > 0:
+ if i > 20:
break
- ic_low = ic
+ Da_low = Da
- print(ic)
+ print(Da)
num_subplots_target = 8
@@ -157,7 +195,7 @@ if __name__ == "__main__":
P = mesh.get_plot_matrix(y_final)
temp_min = 0.0
- temp_max = 0.5
+ temp_max = 5
# Determine subplot grid
ncols = 4
@@ -183,7 +221,7 @@ if __name__ == "__main__":
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_title('T = 0.00 (Initial)')
ax_current.set_xticks([])
ax_current.set_yticks([])
@@ -197,7 +235,7 @@ if __name__ == "__main__":
# 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)
+ 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):
ax_current = axes[plots_created] # Next available subplot
@@ -212,7 +250,7 @@ if __name__ == "__main__":
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_title('T = {:.2f}'.format(current_time))
ax_current.set_xticks([])
ax_current.set_yticks([])
plots_created += 1
@@ -227,5 +265,107 @@ if __name__ == "__main__":
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()
+ fig.suptitle('Temperature Distribution Over Time')
+ 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()