Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • main
1 result

Target

Select target project
  • twd1g21/flash-flood
  • jp7g21/flash-flood
2 results
Select Git revision
  • main
1 result
Show changes
Commits on Source (4)
*~
*.pdf
*.bkp
\ No newline at end of file
This diff is collapsed.
from scipy import integrate as sp_int
from math import exp
from matplotlib import pyplot as plt
import numpy as np
import math
Da = 1
L = 1
delta_x = 0.02
T = 1
delta_t = delta_x**2 / (4*Da)
def step(t, y):
dy_dt = []
for i, yi in enumerate(y):
if i == 0:
yl = 0 #y[0] - (y[1] - y[0]) / delta_x
yr = y[1]
#dyi_dt = -yi
elif i == len(y)-1:
yr = 0 #y[-1] + (y[-1] - y[-2]) / delta_x
yl = y[-2]
#dyi_dt = -yi
else:
yr = y[i+1]
yl = y[i-1]
dyi_dt = (yl - 2*yi + yr) / (Da * delta_x**2)
if yi < 15:
dyi_dt += exp(yi)
else:
pass #raise Exception()
dy_dt.append(dyi_dt)
return dy_dt
def power_gen(y):
P = []
for yi in y:
Pi = 0
for j, yij in enumerate(yi):
if j == 0:
pass
elif j == len(y)-1:
pass
else:
Pi += exp(yij)*delta_x
P.append(Pi)
return P
def power_loss(y):
P = [(yi[1]-yi[0] + yi[-2]-yi[-1])/(Da * delta_x) for yi in y]
return P
a_bounds = [0.53784, 0.537841] #[0.513, 0.6]
t = np.array([0.01*i for i in range(101)])
#y0 = [1.75 for i in range(int(L/delta_x)+1)]
x = np.array([delta_x*i for i in range(int(L/delta_x)+1)])
stddev = 0.05
k = 0
while k < 1:
#a = 0.5 * (a_bounds[0] + a_bounds[1])
a = 0.53784 #0.5378404248361586 #0.53130209261235
print(a)
y0 = a * 1/(stddev*math.sqrt(2*math.pi)) * np.exp(-(x-L/2)**2 / (2*stddev))
try:
sln = sp_int.solve_ivp(step, [0, T], y0, t_eval=t, max_step=delta_t,
atol=1e-2, rtol=1e-3)
a_bounds[0] = a
k += 1
except:
a_bounds[1] = a
for i, y in enumerate(sln.y.T):
if (i % 5) == 0:
#if i == 80:
plt.plot(x, y)
#A = 0.758582299525377
#A = 5.4 #5.46935138606105
#plt.plot(x, np.log(2*A**2/Da) - 2*np.log(np.cosh(A*(x-1.04/2))),
# color="red")
plt.figure()
Pg = power_gen(sln.y.T)
Pl = power_loss(sln.y.T)
P = np.array(Pg) - np.array(Pl)
plt.plot(t, Pg)
plt.plot(t, Pl)
plt.plot(t, P)
plt.show()
from scipy import integrate as sp_int
from math import exp
from matplotlib import pyplot as plt
import numpy as np
from numpy.linalg import inv
import math
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
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
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 discrete_laplacian_2d_cyl(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 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)
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
# def squared_change_rate(y, M):
# dy_dt = np.matmul(M, y) + np.exp(y)
# return np.dot(dy_dt, dy_dt)
# def sln_change(y, M):
# sln, i = solve(M, y, 10, 0.1)
# change = np.abs(sln[10,:] - y) - np.abs(sln[99,:] - y) - 1e6 * i**2
# return (change**2).sum()
# def get_equillibrium(M):
# ymin = minimize(lambda y: sln_change(y, M),
# 1*np.ones(M.shape[0]), method="Nelder-Mead")
# print(ymin)
# return ymin.x
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("")
ic = 0
Da_low = 0
Da_high = 10
i = 0
while i < 100:
i += 1
Da = (Da_low + Da_high) / 2
M = discrete_laplacian_2d_cyl(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(Da)
num_subplots_target = 8
y_final = sln[-1,:]
P = mesh.get_plot_matrix(y_final)
temp_min = 0.0
temp_max = 5
# 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.00 (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 - 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
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}'.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')
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()
......@@ -396,17 +396,11 @@ def main():
}
plots = [
#plot_defaults | { "prefix": "fig/baseline" },
#plot_defaults | { "prefix": "alpha-0.02", "alpha_val": 0.02 },
#plot_defaults | { "prefix": "f-0.1", "f_val": 0.1 },
#plot_defaults | { "prefix": "dx_0.25,dt_0.25", "delta_x": 0.25, "delta_t": 0.25 },
#plot_defaults | { "prefix": "fig/a-5", "a_val": 5 },
#plot_defaults | { "prefix": "fig/a-3", "a_val": 3, "L": 700, "T": 400 },
#plot_defaults | { "prefix": "fig/a-10", "a_val": 10 },
#plot_defaults | { "prefix": "fig/a-1", "a_val": 1 },
#plot_defaults | { "prefix": "fig/slope-0.1", "A_init": gen_A_init(20,1,100,0.1) },
#plot_defaults | { "prefix": "fig/A-pulse", "A_init": 1 + 19*1000/(1*1000+(x_sym-100)**2) },
plot_defaults | { "prefix": "fig/A-pulse-short", "A_init": 1 + 19*10/(1*10+(x_sym-100)**2) }
plot_defaults | { "prefix": "fig/baseline" },
plot_defaults | { "prefix": "fig/alpha-0.02", "alpha_val": 0.02 },
plot_defaults | { "prefix": "fig/a-10", "a_val": 10 },
plot_defaults | { "prefix": "fig/a-1", "a_val": 1 },
plot_defaults | { "prefix": "fig/A-pulse", "A_init": 1 + 19*1000/(1*1000+(x_sym-100)**2) },
]
plt.rcParams.update({'font.size': 22})
......