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

Dunno

parent d86745d1
Branches
No related tags found
No related merge requests found
...@@ -6,6 +6,7 @@ from numpy.linalg import inv ...@@ -6,6 +6,7 @@ from numpy.linalg import inv
import math import math
from scipy.integrate import solve_ivp from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from scipy.optimize import minimize
def sum_to(k): def sum_to(k):
return k*(k+1)/2 return k*(k+1)/2
...@@ -102,8 +103,23 @@ def solve(M, y0, T, delta_t): ...@@ -102,8 +103,23 @@ def solve(M, y0, T, delta_t):
return sln, i+1 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
if __name__ == "__main__": if __name__ == "__main__":
mesh = Mesh2D(50, 50, 40) mesh = Mesh2D(50, 50, 60)
# 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="")
...@@ -112,7 +128,7 @@ if __name__ == "__main__": ...@@ -112,7 +128,7 @@ if __name__ == "__main__":
M = construct_matrix(mesh) / Da M = construct_matrix(mesh) / Da
ic_low = 0 #3.4 ic_low = 0 #3.4
ic_high = 10 #3.4 ic_high = 0 #3.4
i = 0 i = 0
while i < 100: while i < 100:
...@@ -125,10 +141,11 @@ if __name__ == "__main__": ...@@ -125,10 +141,11 @@ if __name__ == "__main__":
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(ic, N_t_in_sln)
if N_t_in_sln < round(simulation_time / time_step): if N_t_in_sln < round(simulation_time / time_step):
ic_high = ic ic_high = ic
else: else:
if i > 50: if i > 0:
break break
ic_low = ic ic_low = ic
...@@ -140,7 +157,7 @@ if __name__ == "__main__": ...@@ -140,7 +157,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 = 6 temp_max = 0.5
# Determine subplot grid # Determine subplot grid
ncols = 4 ncols = 4
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment