请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 369|回复: 12

[源代码] python求解势函数和流函数的简单办法

[复制链接]

新浪微博达人勋

发表于 2024-11-26 21:46:16 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
参考沙颖凯 2014-3-2的帖子 https://bbs.06climate.com/forum. ... =%CA%C6%BA%AF%CA%FD
lz把matlab代码转化成了经常使用的python,发现大佬用的差分格式是五点的空间差分,算法是松弛迭代法,求解泊松方程得到得势函数,在求解时,将边值设置成0视为合理。
虽然改进后得迭代法已经比原来快了一倍,但还是很慢啊……有木有
在等待了N个小时之后,lz意识到泊松方程是线性的,且我们用的差分格式决定了我们的线性系数矩阵是稀疏的,lz尝试了用python线性求解器,非常快速地解决了这个问题。
这里省略了散度场的求解,将散度场作为输入量,求解线性方程直接得到势函数
同理,在求解流函数时,只要把散度场替换成涡度场,求解线性方程也能直接得到流函数
实现了流函数势函数自由~
————————————————————————————————————
欢迎交流讨论~~~

def dx_atmos(longitude, latitude):
    R = 6378137.0                                  # 赤道半径,单位为米
    dtheta = np.gradient(longitude,axis=1)                    # 计算经度的梯度
    dtheta=np.array(dtheta)
    dx = dtheta  * np.cos(np.radians(latitude)) * (np.pi / 180) * R  # 转换为米
    return dx
def dy_atmos(latitude):
    R = 6356752.3                          # 经圈半径,单位为米
    dtheta = np.gradient(latitude,axis=0)             # 计算纬度的梯度
    dtheta = np.array(dtheta)
    dy = dtheta * (np.pi / 180) * R           # 转换为米
    return dy

def grad_atmos(longitude, latitude, H):
    dx = dx_atmos(longitude, latitude)        # 计算经向网格间距
    dy = dy_atmos(latitude)                    # 计算纬向网格间距

    # 计算场 H 的梯度
    grad_y, grad_x = np.gradient(H)           # 计算 H 在经度和纬度方向的梯度

    # 归一化梯度,计算 H 在经度和纬度方向的实际梯度
    grad_x = grad_x / dx                       # 将经度方向的梯度进行归一化
    grad_y = grad_y / dy                       # 将纬度方向的梯度进行归一化

    return grad_x, grad_y

def index(i, j, ny):
    return i * ny + j

def chi_potential(latitude, longitude, divh):
    """
    ps:为了让结果更可靠
    计算速度势和辐散风分量。
    返回:
    chi : np.ndarray
        速度势数组。
    Uchi : np.ndarray
        辐散风 U 分量。
    Vchi : np.ndarray
        辐散风 V 分量。
    """   
    M  = latitude.shape[0]                      # 获取经度和纬度的维度
    N  = latitude.shape[1]
    #divh = divh_atmos(longitude, latitude, u, v)  # 计算水平散度

    dx2 = dx_atmos(longitude, latitude) ** 2   # 计算纬向梯度的平方
    dy2 = dy_atmos(latitude) ** 2               # 计算经向梯度的平方
    #print(dx2.shape)
    #print(f"dy2:{dy2.mean()}")
    #print(f"div:{divh.mean()}")
    #divh*=1e7

    # 2. 构建离散化矩阵
    n_interior = (M - 2) * (N - 2)  # 内部网格点数
    A = sp.lil_matrix((n_interior, n_interior))  # 创建稀疏矩阵
    #print(f"A.shape:{A.shape}")
    # 填充矩阵 A
    for i in range(1,M - 1):
        for j in range(1, N - 1):
            idx = index(i - 1 , j - 1, N - 2)
            #print(f"i:{i} j:{j}")
            A[idx, idx] = 2*((1/dx2[i,j])+(1/dy2[i,j]))  # 对角线
            if i>1:
                A[idx, index(i - 1 - 1, j - 1, N-2)] = -1/dy2[i-1,j]   # 上侧点
            if i<M-2:
                A[idx, index(i + 1 - 1, j - 1, N-2)] = -1/dy2[i+1,j]   # 下侧点
            if j>1:
                A[idx, index(i - 1, j - 1 - 1, N-2)] = -1/dx2[i,j-1]   # 左侧点
            if j<N-2:
                A[idx, index(i - 1, j + 1 - 1, N-2)] = -1/dx2[i,j+1]   # 右侧点

    # 3. 定义右边项
    b = np.zeros(n_interior)

    # 这里定义一个示例源项
    for i in range(1, M - 1):
        for j in range(1, N - 1):
            b[index(i - 1, j - 1 , N - 2)] = divh[i,j]
    # 4. 使用 pyAMG 求解线性系统
    ml = pyamg.ruge_stuben_solver(A)  # 创建多重网格求解器
    chi_interior = ml.solve(b, tol=1e-10)  # 求解

    # 结果转换为二维数组格式用于可视化
    chi = np.zeros((M, N))  # 初始化整个解向量
    chi[1:-1, 1:-1] = chi_interior.reshape((M - 2, N - 2))  # 将内部解放入完整解中

    #print(f"chi:{chi.mean()}")
    #print(f"chi[10,10]:{chi[10,10]}")

    # 计算辐散风分量
    DchiDx, DchiDy = grad_atmos(longitude, latitude, chi)
    Uchi = -DchiDx  # 辐散风的 U 分量
    Vchi = -DchiDy  # 辐散风的 V 分量

    # 将边界值设为 NaN
    chi[0, :] = np.nan; Uchi[0, :] = np.nan; Vchi[0, :] = np.nan
    chi[-1, :] = np.nan; Uchi[-1, :] = np.nan; Vchi[-1, :] = np.nan
    chi[:, 0] = np.nan; Uchi[:, 0] = np.nan; Vchi[:, 0] = np.nan
    chi[:, -1] = np.nan; Uchi[:, -1] = np.nan; Vchi[:, -1] = np.nan

    return chi, Uchi, Vchi



密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2024-11-26 22:08:46 | 显示全部楼层

回帖奖励 +5 金钱

楼主可以参考下windspharm这个库,也很好用
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-11-26 22:18:18 | 显示全部楼层
SunJiaming 发表于 2024-11-26 22:08
楼主可以参考下windspharm这个库,也很好用

windspharm真的有人装上了吗,听说能满足各种气象上计算需求,就是安装了两次都卡在import某个包上了(叹气
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-11-27 00:45:42 | 显示全部楼层
windspharm+1哈哈
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-11-27 08:03:09 | 显示全部楼层
木叶村洗剪吹 发表于 2024-11-26 22:18
windspharm真的有人装上了吗,听说能满足各种气象上计算需求,就是安装了两次都卡在import某个包上了(叹 ...

我是在conda里建了个虚拟环境单独装的,这个库貌似对python的版本有要求
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-11-27 10:29:39 | 显示全部楼层
木叶村洗剪吹 发表于 2024-11-26 22:18
windspharm真的有人装上了吗,听说能满足各种气象上计算需求,就是安装了两次都卡在import某个包上了(叹 ...

windspharm对python和numpy的版本有要求,可以参考这篇帖子下的回复http://bbs.06climate.com/forum.p ... ighlight=windspharm
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-11-27 17:07:52 | 显示全部楼层
GXYYYYY 发表于 2024-11-27 10:29
windspharm对python和numpy的版本有要求,可以参考这篇帖子下的回复http://bbs.06climate.com/for ...

十分感谢!!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-12-2 11:16:44 | 显示全部楼层

如果其他方法都没解决,记得把python版本降到3.9及以下,windspharm没法在3.10以后版本上运行,就是装上了也没用
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-12-2 21:04:49 | 显示全部楼层
本帖最后由 木叶村洗剪吹 于 2024-12-2 21:51 编辑
linktv 发表于 2024-12-2 11:16
如果其他方法都没解决,记得把python版本降到3.9及以下,windspharm没法在3.10以后版本上运行,就是装上 ...

哈哈哈在上面的帖子里眼熟大佬了(!谢谢!已经装上了!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 7 天前 | 显示全部楼层
村里也是出息了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表