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

Added fire code

parent ac9054db
No related branches found
No related tags found
No related merge requests found
*~
*.pdf
*.bkp
\ No newline at end of file
This diff is collapsed.
fire.py 0 → 100644
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
import math
class Mesh2D:
def __init__(self, R, h, delta_r, delta_h):
......@@ -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})
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment