Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
F
flash-flood
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Deploy
Releases
Package registry
Model registry
Operate
Terraform modules
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
GitLab community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
jp7g21
flash-flood
Compare revisions
main to main
Compare revisions
Changes are shown as if the
source
revision was being merged into the
target
revision.
Learn more about comparing revisions.
Source
jp7g21/flash-flood
Select target project
No results found
main
Select Git revision
Branches
main
1 result
Swap
Target
twd1g21/flash-flood
Select target project
twd1g21/flash-flood
jp7g21/flash-flood
2 results
main
Select Git revision
Branches
main
1 result
Show changes
Only incoming changes from source
Include changes to target since source was created
Compare
Commits on Source (4)
Added fire code
· cfb86c86
Joe Pater
authored
1 month ago
cfb86c86
2D case working
· d86745d1
Joe Pater
authored
1 month ago
d86745d1
Dunno
· cd036439
Joe Pater
authored
1 month ago
cd036439
Incorporated 1D case into 2D code
· 42ec5bc6
Joe Pater
authored
1 month ago
42ec5bc6
Expand all
Show whitespace changes
Inline
Side-by-side
Showing
5 changed files
.gitignore
+2
-0
2 additions, 0 deletions
.gitignore
diagrams.drawio
+171
-99
171 additions, 99 deletions
diagrams.drawio
fire.py
+98
-0
98 additions, 0 deletions
fire.py
fire2d.py
+371
-0
371 additions, 0 deletions
fire2d.py
flood.py
+5
-11
5 additions, 11 deletions
flood.py
with
647 additions
and
110 deletions
.gitignore
View file @
42ec5bc6
*~
*.pdf
*.bkp
\ No newline at end of file
This diff is collapsed.
Click to expand it.
diagrams.drawio
View file @
42ec5bc6
This diff is collapsed.
Click to expand it.
fire.py
0 → 100644
View file @
42ec5bc6
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
()
This diff is collapsed.
Click to expand it.
fire2d.py
0 → 100644
View file @
42ec5bc6
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
()
This diff is collapsed.
Click to expand it.
flood.py
View file @
42ec5bc6
...
...
@@ -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
})
...
...
This diff is collapsed.
Click to expand it.