{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# 二维声波问题\n", "\n", "[![下载Notebook](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/resource/_static/logo_notebook.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/r2.5.0/mindflow/zh_cn/cfd_solver/mindspore_acoustic.ipynb) [![下载样例代码](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/resource/_static/logo_download_code.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/r2.5.0/mindflow/zh_cn/cfd_solver/mindspore_acoustic.py) [![查看源文件](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/resource/_static/logo_source.svg)](https://gitee.com/mindspore/docs/blob/r2.5.0/docs/mindflow/docs/source_zh_cn/cfd_solver/acoustic.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 概述\n", "\n", "声波方程求解是医疗超声、地质勘探等领域中的核心技术，大规模声波方程求解面临算力和存储的挑战。声波方程求解器一般采用频域求解算法和时域求解算法，时域求解算法的代表是时域有限差分法 (TDFD)，频域求解算法包括频域有限差分法 (FDFD)、有限元法 (FEM) 和 CBS (Convergent Born series) 迭代法。CBS 方法由于不引入频散误差，且求解的内存需求低，因此受到工程和学术界的广泛关注。尤其是 [Osnabrugge et al. (2016)](https://linkinghub.elsevier.com/retrieve/pii/S0021999116302595) 解决了该方法的收敛性问题，使得 CBS 方法的应用具有更广阔的前景。\n", "\n", "本案例将演示如何调用 MindFlow 提供的 CBS API 实现二维声波方程的求解。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 问题描述\n", "\n", "声波方程求解中，波速场和震源信息是输入参数，求解输出的是时空分布的波场。\n", "\n", "二维声波方程的表达式如下\n", "\n", "| 时域表达式 | 频域表达式 |\n", "| ----------------------------------------------------- | ------------------------------------------------- |\n", "| $\\frac{\\partial^2u}{\\partial t^2} - c^2 \\Delta u = f$ | $-\\omega^2 \\hat{u} - c^2 \\Delta\\hat{u} = \\hat{f}$ |\n", "\n", "其中\n", "\n", "- $u(\\mathbf{x},t) \\;\\; [L]$ 变形位移 (压强除以密度)，标量\n", "- $c(\\mathbf{x}) \\;\\; [L/T]$ 波速，标量\n", "- $f(\\mathbf{x},t) \\;\\; [L/T^2]$ 震源激励 (体积分布力)，标量\n", "\n", "实际求解中，为了降低参数维度，一般先将参数无量纲化，然后针对无量纲方程和参数进行求解，最后恢复解的量纲。选取 $\\omega$、$\\hat{f}$ 和 $d$（网格间距，本案例要求网格在各方向间距相等）对频域方程做无量纲化，可得频域无量纲方程：\n", "\n", "$$\n", "u^* + c^{*2} \\tilde{\\Delta} + f^* = 0\n", "$$\n", "\n", "其中\n", "\n", "- $u^* = \\hat{u} \\omega^2 / \\hat{f}$ 为无量纲变形位移\n", "- $c^* = c / (\\omega d)$ 为无量纲波速\n", "- $\\tilde{\\Delta}$ 为归一化 Laplace 算子，即网格间距均为 1 时的 Laplace 算子\n", "- $f^*$ 为标记震源位置的 mask，即在震源作用点为 1，其余位置为 0\n", "\n", "本案例中 `src` 包可以在 [src](https://gitee.com/mindspore/mindscience/tree/r0.7/MindFlow/applications/cfd/acoustic/src) 下载。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "import numpy as np\n", "import scipy\n", "import pandas as pd\n", "import mindspore as ms\n", "from mindspore import Tensor\n", "\n", "from mindflow.utils import load_yaml_config\n", "\n", "from cbs.cbs import CBS\n", "from src import visual\n", "from solve_acoustic import solve_cbs" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 定义输入参数和输出采样方式\n", "\n", "本案例所需的输入为有量纲 2D 速度场、震源位置列表、震源波形，在文件 [config.yaml](https://gitee.com/mindspore/mindscience/blob/r0.7/MindFlow/applications/cfd/acoustic/config.yaml) 中指定输入文件名。为了方便用户直接验证，本案例在本[链接](https://download-mindspore.osinfra.cn/mindscience/mindflow/dataset/applications/cfd/acoustic)中提供了预置的输入数据，请下载所需要的数据集，并保存在 `./dataset` 目录下。数据集包括速度场 `velocity.npy`、震源位置列表 `srclocs.csv`、震源波形 `srcwaves.csv`。用户可仿照输入文件格式自行修改输入参数。\n", "\n", "输出为时空分布的波场，为了明确输出如何在时间和频率上采样，需在 [config.yaml](https://gitee.com/mindspore/mindscience/blob/r0.7/MindFlow/applications/cfd/acoustic/config.yaml) 文件中指定 `dt`, `nt` 等参数。\n", "\n", "由于输入的震源波形在时间上的采样率可能与输出所要求的不一致，因此需对其进行插值。" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "ms.set_context(device_target='Ascend', device_id=0, mode=ms.GRAPH_MODE)\n", "\n", "config = load_yaml_config('config.yaml')\n", "\n", "data_config = config['data']\n", "solve_config = config['solve']\n", "summary_config = config['summary']\n", "\n", "# read time & frequency points\n", "dt = solve_config['dt']\n", "nt = solve_config['nt']\n", "ts = np.arange(nt) * dt\n", "omegas_all = np.fft.rfftfreq(nt) * (2 * np.pi / dt)\n", "\n", "# read source locations\n", "df = pd.read_csv(os.path.join(data_config['root_dir'], data_config['source_locations']), index_col=0)\n", "slocs = df[['y', 'x']].values # shape (ns, 2)\n", "\n", "# read & interp source wave\n", "df = pd.read_csv(os.path.join(data_config['root_dir'], data_config['source_wave']))\n", "inter_func = scipy.interpolate.interp1d(df.t, df.f, bounds_error=False, fill_value=0)\n", "src_waves = inter_func(ts) # shape (nt)\n", "src_amplitudes = np.fft.rfft(src_waves) # shape (nt//2+1)\n", "\n", "# read velocity array\n", "velo = np.load(os.path.join(data_config['root_dir'], data_config['velocity_field']))\n", "nz, nx = velo.shape\n", "dx = data_config['velocity_dx']" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 选取待求频点\n", "\n", "确定了输出采样方式即确定了所有待求频点。但为了减少计算量，也可以只选择部分频点进行求解，其余频点通过插值获得。具体的频点降采样方式由 [config.yaml](https://gitee.com/mindspore/mindscience/blob/r0.7/MindFlow/applications/cfd/acoustic/config.yaml) 文件中的 `downsample_mode`, `downsample_rate` 指定。默认为不做降采样，即求解除 $\\omega=0$ 之外的所有频点。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# select omegas\n", "no = len(omegas_all) // solve_config['downsample_rate']\n", "\n", "if solve_config['downsample_mode'] == 'exp':\n", " omegas_sel = np.exp(np.linspace(np.log(omegas_all[1]), np.log(omegas_all[-1]), no))\n", "elif solve_config['downsample_mode'] == 'square':\n", " omegas_sel = np.linspace(omegas_all[1]**.5, omegas_all[-1]**.5, no)**2\n", "else:\n", " omegas_sel = np.linspace(omegas_all[1], omegas_all[-1], no)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 执行仿真\n", "\n", "将相关数组定义为 Tensor，调用 `solve_cbs()`，在 NPU 执行求解。由于显存限制，求解过程在频点维度分批执行，batch 数量由用户在 `config.yaml` 中指定，不要求整除频点数（允许最后一个 batch 的 size 与其他 batch 不一致）。求解完成后，会保存频域求解结果到文件 `u_star.npy`。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# send to NPU and perform computation\n", "os.makedirs(summary_config['root_dir'], exist_ok=True)\n", "velo = Tensor(velo, dtype=ms.float32, const_arg=True)\n", "cbs = CBS((nz, nx), remove_pml=False)\n", "\n", "ur, ui = solve_cbs(cbs, velo, slocs, omegas_sel, dx=dx, n_batches=solve_config['n_batches']) # shape (ns, no, len(receiver_zs), nx)\n", "\n", "u_star = np.squeeze(ur.numpy() + 1j * ui.numpy()) # shape (ns, no, len(krs), nx)\n", "np.save(os.path.join(summary_config['root_dir'], 'u_star.npy'), u_star)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 仿真结果后处理\n", "\n", "CBS 求解的是无量纲化的频域方程，但下游任务通常希望在时域观察有量纲波场的演化过程，因此最后将求解结果恢复量纲化并转回时域。恢复量纲的方式为 $\\hat{u} = u^* \\hat{f} / \\omega^2$，若在前述的“选取待求频点”步骤中对频点做了降采样，则在此处需沿频率方向插值恢复所有频点的解。然后对有量纲的频域波场 $\\hat{u}$ 做 Fourier 反变换得到时域波场 $u$。将时域波场保存至文件 `u_time.npy`。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# recover dimension and interpolate to full frequency domain\n", "u_star /= omegas_sel.reshape(-1, 1, 1)**2\n", "u_star = scipy.interpolate.interp1d(omegas_sel, u_star, axis=1, kind='cubic', bounds_error=False, fill_value=0)(omegas_all)\n", "u_star *= src_amplitudes.reshape(-1, 1, 1)\n", "\n", "# transform to time domain\n", "u_time = np.fft.irfft(u_star, axis=1)\n", "np.save(os.path.join(summary_config['root_dir'], 'u_time.npy'), u_time)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "最后，读取时域波场并可视化。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# visualize the result\n", "u_time = np.load(os.path.join(summary_config['root_dir'], 'u_time.npy'))\n", "visual.anim(velo.numpy(), u_time, ts, os.path.join(summary_config['root_dir'], 'wave.gif'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![wave.gif](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/docs/mindflow/docs/source_zh_cn/cfd_solver/images/wave.gif)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1" } } }, "nbformat": 4, "nbformat_minor": 2 }