24小时热门版块排行榜    

查看: 1716  |  回复: 9

Ego悄吟

新虫 (初入文坛)

[求助] 基于有限差分的波场模拟,pml边界,震源为单力源已有2人参与

本人最近在做波场模拟的程序设计,在波场分离的时候遇到下图问题,震源附近发生突变,有没有大佬能够指出问题所在?
模型大小500x500x500,vp为3200km/s,vs为vp/1.732,密度为2250kg/m^3,震源位于中心,采用雷克子波,中心频率为50hz。

基于有限差分的波场模拟,pml边界,震源为单力源
水平.jpg


基于有限差分的波场模拟,pml边界,震源为单力源-1
3D.jpg
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

diwushuxue

铁虫 (初入文坛)

【答案】应助回帖

感谢参与,应助指数 +1
网格多大,你这跟波场分离没关系,是你正演没做好。有可能是频散过于严重
2楼2019-04-17 11:39:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Ego悄吟

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by diwushuxue at 2019-04-17 11:39:28
网格多大,你这跟波场分离没关系,是你正演没做好。有可能是频散过于严重

您好,看到您的回帖之后,我回去划细网格了,之前网格大小为10m,计算步长为1.60e-3s,之后网格大小为2.5m,计算步长为0.8e-3s,可是依然出现这个情况。
基于有限差分的波场模拟,pml边界,震源为单力源-2
3D.jpg


基于有限差分的波场模拟,pml边界,震源为单力源-3
水平.jpg

3楼2019-04-18 09:21:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

diwushuxue

铁虫 (初入文坛)

引用回帖:
3楼: Originally posted by Ego悄吟 at 2019-04-18 09:21:33
您好,看到您的回帖之后,我回去划细网格了,之前网格大小为10m,计算步长为1.60e-3s,之后网格大小为2.5m,计算步长为0.8e-3s,可是依然出现这个情况。

3D.jpg

水平.jpg
...

还是正演问题,你空间网格10m,时间采样1ms,主频30Hz。再试一下!
4楼2019-04-18 17:34:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Ego悄吟

新虫 (初入文坛)

引用回帖:
4楼: Originally posted by diwushuxue at 2019-04-18 17:34:15
还是正演问题,你空间网格10m,时间采样1ms,主频30Hz。再试一下!...

您好,根据您给的建议,我回去试了试,结果还是出现之前的问题。
基于有限差分的波场模拟,pml边界,震源为单力源-4
3D.jpg


基于有限差分的波场模拟,pml边界,震源为单力源-5
水平.jpg

5楼2019-04-19 09:57:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

baobiao007

木虫 (职业作家)

中国特色

明显是正演的递推过程没写对。。
我同意叔本华的观点,人们投身艺术和科学领域的强烈愿望之一就是逃离痛苦、残酷和枯燥无味的现实生活,逃离自己飘忽不定的七情六欲的桎梏。--爱因斯坦
6楼2019-04-19 19:21:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Ego悄吟

新虫 (初入文坛)

引用回帖:
6楼: Originally posted by baobiao007 at 2019-04-19 19:21:44
明显是正演的递推过程没写对。。

您好,您的意思是差分出错了么?我检查了我的程序,貌似并没有问题。

发自小木虫IOS客户端
7楼2019-04-20 08:39:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

baobiao007

木虫 (职业作家)

中国特色

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +1
Ego悄吟: 金币+5, ★★★★★最佳答案 2019-04-21 09:30:08
送你段例子代码吧,直接运行即可


% Finite-differences in time domain(FDTD) acoustic wave propagation in 2D
% medium

close all;
% Output every ... time steps
IT_DISPLAY = 10;

%% MODEL
% Model dimensions
nx = 401;
nz = 401;
dx = 10;    % [m]
dz = 10;    % [m]

% Elastic parameters
vp = 3300.0 * ones(nz, nx);         % velocity of compressional waves, [m/s]

%% TIME STEPPING
t_total = .55;                      % [sec] recording duration
dt = 0.7 * min(dx,dz)/max(vp();   % min grid space / max velocity
nt = round(t_total/dt);             % number of time steps
t = [0:nt]*dt;

CFL = max(vp() * dt / min(dx,dz);
%% SOURCE
f0 = 10.0;                          % dominant frequency of the wavelet
t0 = 1.20 / f0;                     % half Ricker wavelet excitation time
factor = 1e10;                      % amplitude coefficient
angle_force = 90.0;                 % spatial orientation of source
                                    % (not relevant for acoustic case)
jsrc = round(nz/2);                 % source location along OZ
isrc = round(nx/2);                 % source location along OX

a = pi*pi*f0*f0;
dt2 = dt^2;
source_term = factor * exp(-a*(t-t0).^2);                            % Gaussian

force_x = sin(angle_force * pi / 180) * source_term * dt2 / (dx * dz);

min_wavelengh = 0.5*min(vp(vp>330))/f0;     % shortest wavelength bounded by velocity in the air


%% SUMMARY
fprintf('#################################################\n');
fprintf('2D acoustic FDTD wave propagation in isotripic \nmedium in displacement formulation\n');
fprintf('#################################################\n');
fprintf('Model:\n\t%d x %d\tgrid nz x nx\n\t%.1e x %.1e\t[m] dz x dx\n',nz, nx, dz,dx);
fprintf('\t%.1e x %.1e\t[m] model size\n',nx*dx, nz*dz);
fprintf('\t%.1e...%.1e\t[m/s] vp\n', min(vp(), max(vp());
fprintf('Time:\n\t%.1e\t[sec] total\n\t%.1e\tdt\n\t%d\ttime steps\n',t_total,dt,nt);
fprintf('Source:\n\t%.1e\t[Hz] dominant frequency\n\t%.1f\t[sec] index time\n',f0,t0);
fprintf('Other:\n\t%.1f\tCFL number\n', CFL);
fprintf('\t%.2f\t[m] shortest wavelength\n\t%d, %d\t points-per-wavelength OX, OZ\n', min_wavelengh, floor(min_wavelengh/dx), floor(min_wavelengh/dz));
fprintf('#################################################\n');

%% ALLOCATE MEMORY FOR WAVEFIELD
p3 = zeros(nz+2,nx+2);            % Wavefields at t
p2 = zeros(nz+2,nx+2);            % Wavefields at t-1
p1 = zeros(nz+2,nx+2);            % Wavefields at t-2
% Coefficients for derivatives
co_dxx = 1/dx^2;
co_dzz = 1/dz^2;

%% Loop over TIME
tic;
for it = 1:nt
    p3 = zeros(size(p2));
    % Second-order derivatives
    dp_dxx = co_dxx * (p2(2:end-1,1:end-2) - 2*p2(2:end-1,2:end-1) + p2(2:end-1,3:end));
    dp_dzz = co_dzz * (p2(1:end-2,2:end-1) - 2*p2(2:end-1,2:end-1) + p2(3:end,2:end-1));
    % U(t) = 2*U(t-1) - U(t-2) + G dt2/rho;
    p3(2:end-1,2:end-1) = 2.0*p2(2:end-1,2:end-1) - p1(2:end-1,2:end-1) + (vp.^2).*(dp_dxx + dp_dzz).*dt2;
    % Add source term
    p3(jsrc, isrc) = p3(jsrc, isrc) + force_x(it);
    % Exchange data between t-2 (1), t-1 (2) and t (3) and apply ABS
    p1 = p2 .* weights;
    p2 = p3 .* weights;
    % Output
    if mod(it,IT_DISPLAY) == 0
        fprintf('Time step: %d \t %.4f s\n',it, single(t(it)));
        imagesc(p3); colorbar;
        axis equal tight; colormap jet;
        title(['Step = ',num2str(it),'/',num2str(nt),', Time: ',sprintf('%.4f',t(it)),' sec']);
        drawnow;
    end
end
toc; disp('End');
我同意叔本华的观点,人们投身艺术和科学领域的强烈愿望之一就是逃离痛苦、残酷和枯燥无味的现实生活,逃离自己飘忽不定的七情六欲的桎梏。--爱因斯坦
8楼2019-04-20 10:36:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

baobiao007

木虫 (职业作家)

中国特色

【答案】应助回帖

表情的地方  是冒号 :
我同意叔本华的观点,人们投身艺术和科学领域的强烈愿望之一就是逃离痛苦、残酷和枯燥无味的现实生活,逃离自己飘忽不定的七情六欲的桎梏。--爱因斯坦
9楼2019-04-20 10:38:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Ego悄吟

新虫 (初入文坛)

引用回帖:
8楼: Originally posted by baobiao007 at 2019-04-20 10:36:59
送你段例子代码吧,直接运行即可


% Finite-differences in time domain(FDTD) acoustic wave propagation in 2D
% medium

close all;
% Output every ... time steps
IT_DISPLAY = 10;

%% MODEL
% M ...

您好,谢谢您的代码,我的模型的问题确实出现在差分上。

发自小木虫IOS客户端
10楼2019-04-21 09:28:27
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 Ego悄吟 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 化学B口多少分能上会呀 +4 WOWO159357 2024-05-22 6/300 2024-05-22 18:49 by ckck9999
[教师之家] 博士去高校就是为了有寒暑假吗? +23 wenwen0825 2024-05-16 29/1450 2024-05-22 15:57 by iwdw2012
[考博] 想被211以上高校课题组接收 +11 风起沧澜 2024-05-16 13/650 2024-05-22 14:37 by keyaner23
[硕博家园] 人生 +7 暮色恋伊人 2024-05-22 7/350 2024-05-22 13:52 by 白菜1992
[硕博家园] 博士复试,申请成绩复核,有机会翻盘吗? +15 长海二声笑 2024-05-21 22/1100 2024-05-22 12:44 by 带甲三千
[基金申请] 国社科申请书上传有误,学校已提交到省里,省里还未审核,还能退回修改嘛? 100+3 远山晴岚 2024-05-19 7/350 2024-05-22 12:26 by holypower
[论文投稿] 关于通讯作者 5+3 irikiar 2024-05-21 3/150 2024-05-22 11:21 by bobvan
[基金申请] 国自然等 80+4 胖虎 2024-05-21 12/600 2024-05-22 09:47 by nono2009
[教师之家] 白天不懂夜的黑,90后青椒与60后老板 +6 zylfront 2024-05-18 7/350 2024-05-22 09:28 by songwz
[硕博家园] 耐高温垫片求购 +8 Sexyflea 2024-05-16 11/550 2024-05-22 08:12 by tan151646
[硕博家园] 博三一直没文章怎么办 +27 133456 2024-05-17 45/2250 2024-05-22 06:56 by dong5391
[找工作] 浙江的高校现在门槛都这么高吗 +12 dadqweq_qw 2024-05-16 12/600 2024-05-21 22:30 by foolishmani
[考博] 本科出身不好是不是会被直接刷呀 +5 未来富婆蛙 2024-05-20 5/250 2024-05-21 15:52 by chemdl
[基金申请] 听说面青地E09已经送了么? +6 叉烧吃叉烧 2024-05-21 8/400 2024-05-21 15:42 by 杜大豆豆
[博后之家] 山东大学(青岛)“天然药物生物智造”课题组 招聘“博士后”(年薪20.4-55.6万元) +3 第二种态度 2024-05-18 6/300 2024-05-21 15:37 by 安小樱
[考博] 双非博士还是985科研助理,然后再读博 +6 lxdatj123 2024-05-18 13/650 2024-05-21 08:19 by lxdatj123
[基金申请] 去年申请基金的评审意见是ChatGPT在国内是禁止的,研究方案中有使用ChatGPT不合理 +4 瞬息宇宙 2024-05-19 4/200 2024-05-21 00:37 by dxcharlary
[找工作] 绍兴文理学院怎么样?有没有坑啊 +9 zhaojiang427 2024-05-16 23/1150 2024-05-20 21:44 by 高敖曹
[论文投稿] 求推荐期刊 20+3 好困好困a 2024-05-18 4/200 2024-05-19 11:30 by nono2009
[硕博家园] 五氯化铌怎么溶解啊 +3 南南枝枝 2024-05-17 5/250 2024-05-17 11:37 by ad_fish
信息提示
请填处理意见