-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathTest_Gen_Floor4.py
114 lines (89 loc) · 3.86 KB
/
Test_Gen_Floor4.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
import torch
import os
import numpy as np
import math
import matplotlib.pyplot as plt
from matplotlib import cm
from torch.utils.data import Dataset, DataLoader
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
os.environ['CUDA_VISIBLE_DEVICES'] = '0'
def get_phy_Loss(output1, c, size, loc_x, loc_y, dt, dx, fre):
print(c)
t_max = dt*(size-1)
dy = dx
r1, r2 = c ** 2 * dt ** 2 / dx ** 2, c ** 2 * dt ** 2 / dy ** 2
sum = 0
for n in range(1, int(t_max / dt)):
# 在边界处设置固定边界条件
# 在内部节点上使用五点差分法计算新的波场
output1[(n + 1), 1:-1, 1:-1] = 2 * output1[n, 1:-1, 1:-1] - output1[n - 1, 1:-1, 1:-1] + \
r1[1:-1, 1:-1] * (output1[n, 2:, 1:-1] - 2 * output1[n, 1:-1, 1:-1] + output1[n,:-2, 1:-1]) + \
r2[1:-1, 1:-1] * (output1[n, 1:-1, 2:] - 2 * output1[n, 1:-1, 1:-1] + output1[n,1:-1, :-2])
for idx in range(len(loc_x)):
output1[n + 1, loc_x[idx], loc_y[idx]] = 1500 * np.sin(2 * 3.1415926 * fre * (n + 1) * dt)
output1[n + 1, 0, :] = 0
output1[n + 1, :, 0:1] = output1[n, :, 0:1] - dt * c[:, 0:1] * (output1[n, :, 0:1] - output1[n, :, 1:2]) / dx
output1[n + 1, :, -1:] = output1[n, :, -1:] - dt * c[:, -1:] * (output1[n, :, -1:] - output1[n, :, -2:-1]) / dx
output1[n + 1, -1, :] = output1[n, -1, ] - dt * c[-1, :] * (output1[n, -1:, :] - output1[n, -2:-1, :]) / dx
# output1[:, -1, :] = 0
print(sum)
return torch.unsqueeze(output1, dim=1)
size= 256 + 2
dt = float(1/4096.0)
dx = 1
fre = 25
n = 2
Lx = Ly = 64*n # Length of the 2D domain
loc_x = [37.5*2,37.5*2,37.5*2,37.5*2,37.5*2,37.5*2,37.5*2,37.5*2,37.5*2,
32.5*2,32.5*2,32.5*2,32.5*2,32.5*2,32.5*2,32.5*2,32.5*2,32.5*2,
45*2,45*2,45*2,45*2,45*2,45*2,45*2,45*2,45*2]
loc_y = [4,16,32,48,64,80,96,112,124
,4,16,32,48,64,80,96,112,124
,5,10,15,20,25,106,112,118,124]
np.save('./case/SeaFloor3/loc_x.npy', loc_x)
np.save('./case/SeaFloor3/loc_y.npy', loc_y)
tdx = 1 # 空间步长为0.01米
tdy = 1
c = 1500 * torch.ones((64*n, 64*n)).cuda() # 生成一个波速为45的张量（可根据测试需要调整）
# for i in range(1, c.shape[1] - 1, 1):
# c[i, c.shape[1] - i:] = 0
for i in range(64*n):
if i<32*n:
c[40*n+int(((32*n-i)**2)/((32*n) **2)*(24*n)):,i] = 0
elif i<64*n:
c[40*n + int(((32*n - i) ** 2) / ((32*n) ** 2) * (24*n)):, i] = 0
cmap = cm.get_cmap('jet')
plt.imshow(c.detach().cpu().numpy().squeeze(), cmap=cmap)
plt.colorbar()
plt.savefig('./case/SeaFloor3/speed.png')
plt.show()
#设置时间步
#设置单点（单点声速预测）或多点声源（探测多障碍物）
x1=[[int(50*n),int(37.5*n),int(35*n),int(35*n),int(37.5*n),int(50*n)]]
y1=[[int(6*n),15*n,28*n,38*n,53*n,60*n]]
np.save('./case/SeaFloor3/x1.npy', x1)
np.save('./case/SeaFloor3/y1.npy', y1)
location = torch.ones((len(x1),2)).cuda()
output_s = torch.zeros((size*len(x1), 1, 64*n, 64*n)).cuda()
for i in range(len(x1)):
u = np.zeros((2, 64*n, 64*n))
output = torch.zeros((size, 64*n, 64*n)).cuda()
x = np.arange(0, Lx, tdx)
y = np.arange(0, Ly, tdy)
X, Y = np.meshgrid(x, y)
for idx in range(len(x1[i])):
u[0, x1[i][idx], y1[i][idx]] = 1500 * np.sin(2 * 3.1415926 * fre * 0 * dt)
u[1, x1[i][idx], y1[i][idx]] = 1500 * np.sin(2 * 3.1415926 * fre * 1 * dt)
u = torch.from_numpy(u)
output[0:2,:,:] = u[0:2,:,:]
output=get_phy_Loss(output,c,size,x1[i],y1[i],dt,dx,fre)
output_s[size*i:size*(i+1),:,:,:]=output[:,:,:,:].clone()
cmap = cm.get_cmap('jet')
plt.imshow(output_s[-1].detach().cpu().numpy().squeeze(), cmap=cmap)
plt.colorbar()
plt.show()
#保存文件
torch.save(c,'./case/SeaFloor3/ref_speed.pt')
print(output_s.shape)
torch.save(output_s[:],'./case/SeaFloor3/o_temp.pt')
# torch.save(location,'./case/location.pt')