Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import matplotlib.pyplot as plt
- from math import sin, pi
- import time
- import numba
- @numba.njit
- def q(L, dx, x, y, z, t):
- center = ( (0.2 * L, 0.5 * L, 0.5 * L), (0.5 * L, 0.2 * L, 0.5 * L) )
- size = ( (4 / dx, 0.2 * L, 0.2 * L), (0.2 * L, 4 / dx, 0.2 * L) )
- p = (100, 100)
- tetta = (50, 77)
- fi = (0, 10)
- for i in range(2):
- cx, cy, cz = center[i]
- lx, ly, lz = size[i]
- if (cx - lx <= x <= cx + lx) and (cy - ly <= y <= cy + ly) and (cz - lz <= z <= cz + lz):
- return p[i] * (sin(2 * pi * t / tetta[i] + fi[i]) + 1)
- return 0.0
- @numba.njit
- def D(x, y, z, L):
- """
- Считаем коэффициент дифузии
- Если точка лежит внутри сферы или не лежит
- """
- if 25 * ((x - L / 2) ** 2 + (y - L / 2) ** 2 + (z - L / 2) ** 2) <= L ** 2 :
- return 10**-16 # np.longdouble(10**-16)
- return 1.0 # np.longdouble(1.0)
- def save1DToFile(data, fileName):
- """
- Сохранить в файл одномерный массив:
- - в первой строке размер
- - далее по одному элементу в строке
- """
- with open(fileName, "w") as file:
- file.write(str(len(data)) + '\n')
- for x in data:
- file.write(str(x) + '\n')
- def save2DToFile(data, fileName):
- """
- Сохранить в файл двумерный массив:
- - в первой строке размер
- - далее по одному элементу в строке
- """
- with open(fileName, "w") as file:
- file.write("%d %d\n" % (len(data), len(data[0])))
- for row in data:
- for x in row:
- file.write(str(x) + '\n')
- def save3DToFile(data, fileName):
- """
- Сохранить в файл трехмерный массив:
- - в первой строке размер
- - далее по одному элементу в строке
- """
- with open(fileName, "w") as file:
- file.write("%d %d %d\n" % (len(data), len(data[0]), len(data[0][0])))
- for m in data:
- for row in m:
- for x in row:
- file.write(str(x) + '\n')
- @numba.njit
- def solve(dt, dx, Tm, L, cprobe):
- """
- Выполнение расчетов для заданных параметров
- """
- iter_time = int(Tm / dt) + 1
- iter_dim = int(L / dx) + 1
- T = np.empty((2, iter_dim, iter_dim, iter_dim))
- T.fill(0)
- g = dt / dx ** 2
- probe = np.empty((4, iter_time))
- probe.fill(0)
- cpi = [ [ int(_ / dx) for _ in cprobe[i] ] for i in range(4)] # переводим в индексы из координат
- for t in range(1, iter_time):
- i = t & 1
- j = 1 - i
- # для всех точек кроме граней
- for ix in range(1, iter_dim-1):
- for iy in range(1, iter_dim-1):
- for iz in range(1, iter_dim-1):
- x = ix * dx
- y = iy * dx
- z = iz * dx
- T[i][ix][iy][iz] = (T[j][ix][iy][iz] +
- (g * D(x, y, z, L)) * (
- T[j][ix-1][iy][iz] + T[j][ix+1][iy][iz] +
- T[j][ix][iy-1][iz] + T[j][ix][iy+1][iz] +
- T[j][ix][iy][iz-1] + T[j][ix][iy][iz+1] -
- 6 * T[j][ix][iy][iz] ) +
- dt * q(L, dx, x, y, z, t))
- # на гранях
- # x = L
- ix = iter_dim - 1
- for iy in range(1, iter_dim-1):
- for iz in range(1, iter_dim-1):
- T[i][ix][iy][iz] = T[i][ix-1][iy][iz]
- # z = L
- iz = iter_dim - 1
- for iy in range(1, iter_dim-1):
- for ix in range(1, iter_dim-1):
- T[i][ix][iy][iz] = T[i][ix][iy][iz-1]
- # y = 0
- for ix in range(1, iter_dim-1):
- for iz in range(1, iter_dim-1):
- T[i][ix][0][iz] = T[i][ix][1][iz]
- # z = 0
- for ix in range(1, iter_dim-1):
- for iy in range(1, iter_dim-1):
- T[i][ix][iy][0] = T[i][ix][iy][1]
- # пробы
- for k in range(4):
- ix, iy, iz = cpi[k]
- probe[k][t] = T[i][ix][iy][iz]
- return probe, T[ 1 - (iter_time & 1) ]
- def printColorMap(fileName, L, dx, cprobe):
- """
- Цветовая карта в сечении 0.5 * L в конечный момент времени
- """
- z = int(0.5 * L / dx)
- # считываем данные из файла
- with open(fileName) as file:
- X, Y, Z = map(int, file.readline().split())
- data = [ [ 0 for j in range(Y) ] for i in range(X) ]
- for ix in range(X):
- for iy in range(Y):
- for iz in range(Z):
- v = float(file.readline())
- if iz == z:
- data[ix][iy] = v
- # data = data[::-1][1:]
- data = list(zip(*data))[::-1]
- Y -= 1
- # подготавливаем карту
- fig, ax = plt.subplots(figsize=(40, 8))
- plt.imshow(data, cmap=plt.get_cmap('coolwarm'), extent=[0, X, 0, Y])
- ax.set_aspect(2)
- stepX = max(1, int(4 * (X - 1) / L)) * 4
- stepY = max(1, int(2 * (X - 1) / L)) * 2
- ticksForX = [_ for _ in range(0, X, stepX)]
- ticksLabelForX = ["{:1.2f}".format(_ * dx) for _ in range(0, X, stepX)]
- ticksForY = [_ for _ in range(0, X, stepY)]
- ticksLabelForY = ["{:1.2f}".format(_ * dx) for _ in range(0, X, stepY)]
- ax.set_yticks(ticksForY)
- ax.set_yticklabels(ticksLabelForY)
- ax.set_xticks(ticksForX)
- ax.set_xticklabels(ticksLabelForX)
- ax.set_xlabel("ось X")
- ax.set_ylabel("ось Y")
- # пробы в виде точек на цветовой карте
- for i in range(len(cprobe)):
- x, y, z = cprobe[i]
- ix, iy = int(x/dx), int(y/dx)
- ax.plot([ix], [iy], 'o', color='black')
- ax.annotate(f'({i})', xy=(ix, iy), xytext=(ix+1, iy))
- plt.colorbar()
- def loadProbesFromFiles(probeFile):
- """
- Читаем данные с температурами в пробах из файла для заданного случая case
- """
- probes = [ ]
- for i in range(4):
- with open(probeFile % (i)) as file:
- n = int(file.readline())
- probes.append([])
- for _ in file:
- probes[i].append(float(_))
- return probes
- def showProbes(probeFile, Tm, L, dx, dt):
- """
- Отобразить график изменения температуры в пробах
- """
- # считываем данные из файлов
- probes = loadProbesFromFiles(probeFile)
- n = len(probes[0])
- # рисуем графики
- fig, ax = plt.subplots(figsize=(10, 10))
- t = [_ for _ in range(n)]
- for i in range(4):
- plt.plot(t, probes[i], label=f'проба{i}')
- ticksForX = [_ for _ in range(0, n, (n+9)//10)]
- ticksLabelForX = ["{:1.2f}".format(_ * dt) for _ in range(0, n, (n+9)//10)]
- ax.set_xticks(ticksForX)
- ax.set_xticklabels(ticksLabelForX)
- plt.xlabel("Ось времени (c)")
- plt.ylabel("Температура")
- plt.legend()
- plt.show()
- def doFFT(probeFile, Tm, dt):
- """
- Восстановление периодов колебаний с использованием быстрого преобразования Фурье
- """
- from scipy.signal import find_peaks
- # считываем данные из файлов
- probes = loadProbesFromFiles(probeFile)
- # fig, ax = plt.subplots(figsize=(10, 10))
- for i in range(4):
- sig = np.array(probes[i], dtype=float)
- sig_fft = np.fft.rfft(sig) # строим быстрое преобразование Фурье
- freq = np.fft.rfftfreq(sig.size, d=dt)
- peaks, _ = find_peaks(np.abs(sig_fft))
- print(f'проба{i}, периоды:', end=' ')
- for p in peaks:
- print(f'{1 / freq[p] / dt}', end=' ')
- print()
- def calculateAll(params, L, Tm, cprobe):
- """
- Посчитать все случаи заданные в `params` и сохранить данные в файлы
- или загрузить данные из файлов для отображения
- """
- import os.path
- cp = np.empty((len(cprobe), 3))
- for i in range(len(cprobe)):
- for j in range(3):
- cp[i][j] = cprobe[i][j]
- for i in range(len(params)):
- dx, dt = params[i]
- # проверить, если нет файла расчета - то выполнить
- # dirPath = 'sample_data/save_data/'
- heatFieldFile = dirPath + f'heat-field-case{i+1}.txt'
- probeFile = dirPath + f"probe%d-case{i+1}.txt"
- timeFile = dirPath + f"time-case{i+1}.txt"
- if not os.path.isfile(timeFile):
- startTime = time.time()
- probe, T = solve(dt, dx, Tm, L, cp)
- endTime = time.time()
- # сохранение в файл времени расчета
- calcTime = f"dx = {dx}, dt = {dt}, calculation time = {endTime-startTime} с."
- with open(timeFile, "w") as file:
- file.write(calcTime)
- # сохранение в файлы температуры в пробах
- for j in range(len(probe)):
- save1DToFile(probe[j], probeFile % (j))
- # сохранение в файл температурного поля в конечный момент времени
- save3DToFile(T, heatFieldFile)
- else:
- # если файлы расчетов есть, то не выполнять расчет
- # вывод времени расчета для предыдущего случая
- with open(timeFile) as file:
- calcTime = file.readline().strip()
- # вывод времени расчета
- print(calcTime)
- # построение цветовой карты
- printColorMap(heatFieldFile, L, dx, cprobe)
- # построение графиков проб
- showProbes(probeFile, Tm, L, dx, dt)
- # восстановление периодов
- doFFT(probeFile, Tm, dt)
- params = ((4, 0.2), (2, 0.1), (1, 0.05), (0.5, 0.025))
- L = 100
- Tm = 600
- cprobe = ( (0.3 * L, 0.3 * L, 0.5 * L), (0.3 * L, 0.7 * L, 0.5 * L), (0.7 * L, 0.7 * L, 0.5 * L), (0.7 * L, 0.3 * L, 0.5 * L) )
- calculateAll(params, L, Tm, cprobe)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement