柱坐标下多重网格法解泊松方程-python

写了一个柱坐标下多重网格法解泊松方程的code，外边界采用的是第一类边界条件，直接用解析解赋值，内边界需要注意的是在x=0的转轴上面有奇点。

柱坐标下泊松方程形式为

$\frac{\partial^2 u }{\partial x^2}+\frac{\partial^2 u }{\partial z^2}+\frac{1}{x}\frac{\partial u}{\partial x}=b(x,z)$

当x趋向于0时， $lim\frac{1}{x}\frac{\partial u}{\partial x}=\frac{\partial^2 u }{\partial x^2}$

将方程离散化，其中第三项采用中心差分格式。

开始采用了有残差的多重网格形式，一直不收敛，后来索性残差去掉了，收敛的还可以。

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
import time
def coarsen(x_fine, y_fine, phi_fine):# Lagrange type interpolation, two orderx2 = np.linspace(x_fine[0],x_fine[-1],int(x_fine.size/2))y2 = np.linspace(y_fine[0],y_fine[-1],int(y_fine.size/2))f = interpolate.interp2d(x_fine, y_fine, phi_fine, kind='cubic')#{‘linear’, ‘cubic’, ‘quintic’}phi2 = f(x2,y2)return phi2, x2, y2def fine(x_coarse, y_coarse, phi_coarse):# Lagrange type interpolation, two orderx2 = np.linspace(x_coarse[0],x_coarse[-1],int(x_coarse.size*2))y2 = np.linspace(y_coarse[0],y_coarse[-1],int(y_coarse.size*2))f = interpolate.interp2d(x_coarse, y_coarse, phi_coarse, kind='cubic')#{‘linear’, ‘cubic’, ‘quintic’}phi2 = f(x2,y2)return phi2, x2, y2
def Relax2( b,  phi, x, z,leveli):omega = 1.95ite  = 10*3**leveli   # the iteration time, you can set by other conditionk  =  0h_r = x[1]-x[0]h_z = z[1]-z[0]while(k<ite): #控制迭代次print('relax,k and leveli=',k,leveli)for i in range(0,x.size-1):for j in range(0, z.size-1):phi_ji =  np.copy(phi[j,i])dr_2 = h_r*h_r # dr^2dz_2 = h_z*h_zif (i!=0 and j!=0):# the inner boundarya1 = -2/dr_2-2/dz_2a2 = 1/dr_2+1/x[i]/2/h_ra3 = 1/dr_2-1/x[i]/2/h_ra4 = 1/dz_2a5 = 1/dz_2phi[j,i] = (1.-omega)*phi_ji-omega/a1*(a2* (phi[j,i+1])+a3* (phi[j,i-1])+a4* (phi[j+1,i])+a5* (phi[j-1,i])-b[j,i])continueelif (i!=0 and j==0):a1 = -2/dr_2-2/dz_2 #factor to u[j,i]a2 = 1/dr_2+1/x[i]/2/h_r#factor to u[j,i+1]a3 = 1/dr_2-1/x[i]/2/h_r                #factor to u[j,i-1]a4 = 2/dz_2                #factor to u[j+1,i]phi[j,i] = (1.-omega)*phi_ji-omega/a1*(a2* (phi[j,i+1])+a3* (phi[j,i-1])+a4* (phi[j+1,i])-b[j,i])continueelif (i==0 and j!=0):#z axisa1 = -4/dr_2-2/dz_2a2 = 4/dr_2a3=1/dz_2a4=1/dz_2phi[j,i] = (1.-omega)*phi_ji-omega/a1*(a2* (phi[j,i+1])+a3* (phi[j-1,i])+a4* (phi[j+1,i])-b[j,i])continueelse:a2 = 4/dr_2a3 = 2/dz_2a1 = -1*a2-a3phi[j,i] = (1.-omega)*phi_ji-omega/a1*(a2* (phi[j,i+1])+a3* (phi[j+1,i])-b[j,i])k = k+1return phi
def residual(b, phi, x, z): res  =  np.zeros((z.size,x.size),dtype=float)for i in range(0,x.size-1):for j in range(0,z.size-1):dr_2 = h_r*h_r # dr^2dz_2 = h_z*h_zif (i==0 and j==0):res[j,i] =b[0,0]-4*(phi[0,1]-phi[0,0])/dr_2-2*(phi[1,0]-phi[0,0])/dz_2elif (i!=0 and j==0):a1 = -2/dr_2-2/dz_2 #factor to u[j,i]a2 = 1/dr_2+1/x[i]/2/h_r#factor to u[j,i+1]a3 = 1/dr_2-1/x[i]/2/h_r                #factor to u[j,i-1]a4 = 2/dz_2                #factor to u[j+1,i]res[j,i] = b[j,i]-a1*phi[j,i]-a2*phi[j,i+1]-a3*phi[j,i+1]-a4*phi[j+1,i]elif (i==0 and j!=0):#z axisa1 = -4/dr_2-2/dz_2a2 = 4/dr_2a3=1/dz_2a4=1/dz_2res[j,i]= b[j,i]-a1*phi[j,i]-a2*phi[j,i+1]-a3*phi[j+1,i]-a4*phi[j-1,i]else:a1 = -2/dr_2-2/dz_2a2 = 1/dr_2+1/x[i]/2/h_ra3 = 1/dr_2-1/x[i]/2/h_ra4 = 1/dz_2a5 = 1/dz_2res[j,i] = b[j,i]-a1*phi[j,i]-a2*phi[j,i+1]-a3*phi[j,i-1]-a4*phi[j+1,i]-a5*phi[j-1,i]return resdef bij(x,z):b0 = density_nfw(x,z)return b0
def density_nfw(x,z):                #density of the NFW dark matter potentialrho0_nfw =1# 0.00854*1e9#M0/Kpc^3rh = 19.6 #Kpcrho=np.zeros((z.size,x.size))for j in range(0,z.size):for i in range(0,x.size):r = math.sqrt(x[i]*x[i]+z[j]*z[j])rrh = r/rhif rrh<0.02:rrh=0.02rho[j,i] = rho0_nfw/rrh/(1+rrh)/(1+rrh)return rho
def phi_nfw(x,z):       # the analysical resultrho0_nfw=1#assum 4pi g rho=1rh = 19.6phi0=np.zeros((z.size,x.size),dtype=float)for j in range(0,z.size):for i in range(x.size):r = math.sqrt(x[i]*x[i]+z[j]*z[j])rrh = r/rhif rrh<0.02:rrh=0.02phi0[j,i] =  -1*rh* rh*math.log(1+rrh)/rrhreturn phi0def MG2( b,phi, x,z):h_r0  =  x[1]-x[0] h_z0  =  z[1]-z[0]               vh  =  Relax2(b,  phi, x, z,1)total_level=int(math.log(x.size, 2)-1)t = 0while (t <1): # number of V circles, here is 1for i in range(1, total_level): k = iv2h,x2,z2 = coarsen(x,z,vh)b2 = bij(x2,z2)vh = Relax2(b2,  v2h, x2, z2, k)x = x2z = z2for i in range(1,total_level): v2h,x2,z2 = fine(x,z,vh)b2=bij(x2,z2)vh0 = Relax2(b2,  v2h, x2, z2, k)vh = vh0x = x2z = z2k=k-1t = t+1print('Vtime=',t)return vh'''
def MG2( b,phi, x,z):h_r0  =  x[1]-x[0] h_z0  =  z[1]-z[0]               vh0  =  Relax2(b,  phi, x, z,1)rh  =  residual(b, phi, x, z)# residualeh0 = np.zeros(rh.shape,dtype = float)eh =  Relax2(rh, eh0,  x,z,1)vh = vh0+ehrh  =  residual(b, vh,  x, z)t = 0while ((np.max(rh) > 1e-2) and (t <1)): for i in range(1,5): k = i
#            r2h,x2,z2 =  coarsen(x,z,rh)   # residual to fine gridv2h,x2,z2 = coarsen(x,z,vh)b2=bij(x2,z2)vh0 = Relax2(b2,  v2h, x2, z2, k)rh = residual(b2, vh0, x2, z2)eh =  Relax2(rh, eh0,  x,y,5)vh = vh0+eheh0 = np.zeros(rh.shape,dtype = float)e2h0 = np.zeros(r2h.shape,dtype = float)e2h  =  Relax2(r2h, e2h0,x2,z2,k)phi = v2h+e2hb = bij(x2,z2)vh =  Relax2(b, phi, x2,z2,k)rh = residual(b, vh, x2,z2)x = x2z = z2for i in range(1,5): k = 6-ir2h,x2,z2 =  fine(x,z,rh)   # residual to fine gride2h0 = np.zeros(r2h.shape,dtype = float)# initiate the erre2h  =  Relax2(r2h, e2h0, x2,z2,k)     #calculate errv2h,x2,z2= fine(x,z,vh)                  #to fine gridphi = v2h+e2h                        # update phib = bij(x2,z2)                         # calculate b on this levelvh =  Relax2(b, phi, x2,z2,k)          # calculate  vh again with updated phirh =  residual(b, vh, x2,z2)            #calculate residualx = x2z = z2t = t+1print('Vtime=',t)return vh
''''''
def MG2(b,  phi, x, z):h_r  =  x[1]-x[0]y_r  = z[1]-z[0]vh  =  Relax2(b,  phi, x, z)rh  =  residual(b,  phi, x, z)# residuali = 0while (i<10): #print(phi)r2h,x2 ,z2=  coarsen(x,z,rh)   # residual to coarse gridh_r2=x2[1]-x2[0]h_z2=z2[1]-z2[0]e2h0 = np.zeros((z2.size,x2.size),dtype = float)e2h  =  Relax2(r2h, e2h0,x2,z2)print(e2h.shape)eh,x, z = fine(x2,z2,e2h)print(vh.shape,eh.shape)v1h = vh+ehh_r = x[2]-x[1]h_z = z[1]-z[0]phi = Relax2(b, v1h, x,z )vh = v1hrh = residual(b, vh, x,z)i = i+1print(i)return phi
'''
def initiate():time_start = time.time()n_x_grid =128n_z_grid=128xs = 0.xe =200 # in kpczs=0.ze=200x=np.linspace(xs,xe,n_x_grid)z=np.linspace(zs,ze,n_z_grid)phi  =  np.ones((z.size,x.size), dtype=float)phi2 = phi_nfw(x,z)phi[-1,:] =  (phi2[-1,:])phi[:,-1] =  (phi2[:,-1])b = density_nfw(x,z)h_r=x[1]-x[0]h_z=z[1]-z[0]'''vh0  =  Relax2(b,  phi, x, z,5)rh  =  residual(b, phi, x, z)# residualeh0 = np.zeros(rh.shape,dtype = float)eh =0#  Relax2(rh, eh0,  x,z,5)vh = vh0+eh#rh  =  residual(b, vh,  x,z)'''result =  MG2( b,  phi, x, z) # solve the equationgrav=6.67e-8 #cgsm0=1.989e33kpc=3.086e21rho_halo0=0.00854*1e9*m0/kpc#a^2里面有个kpc的平方potential = 4*math.pi*grav*rho_halo0*resultphi_error = 4*math.pi*grav*rho_halo0*(phi2-result)time_end = time.time() #print(result.shape)time_c= time_end - time_start   #运行所花时间print('time cost', time_c, 's')return potential,phi_error, x, zdef plot(x, z, result,name):print(name)plt.figure(1)plt.contourf(x,z,result ,100,cmap='seismic')plt.colorbar()plt.title(name)plt.xlabel('x(kpc)')plt.ylabel('z(kpc)')plt.savefig(name+'.png')plt.show()plt.close()return 0def plot1d(x,phi,name):print(name)plt.figure(1)plt.plot(x,phi)plt.title(name)plt.xlabel('x(kpc)')plt.ylabel('phi')plt.savefig(name+'.png')plt.show()plt.close()return 0
potent, err, x, z = initiate()
plot(x,z, potent,'Numerical results(erg)')plot(x,z,err,'error')
plot1d(x,err[0,:],'err1d')

柱坐标下多重网格法解泊松方程-python相关推荐

高等数学学习笔记——第七十九讲——柱坐标下三重积分的计算
一. 问题的引入--在直角坐标系下,某些三重积分的计算并不方便,如宇宙飞船体积的计算图略二. 空间上点的柱坐标表示 1. 柱坐标 2. 直角坐标与柱坐标的关系 3. 柱坐标系中的坐标面三. 柱坐 ...
多重网格法解泊松方程（两步法）
这里是一个一维泊松方程,采用两步法求解,教程见 http://bender.astro.sunysb.edu/classes/numerical_methods/lectures/elliptic-m ...
电磁波在柱坐标和球坐标下的本征解
一.背景课程和习题中,我们通常接触的都是平面电磁波,但现实生活中却常常碰到柱面波和球面波,如通电导线辐射场.手机信号等.而且工程上,平面波也可以按柱面波和球面波展开. 那么电磁波在柱坐标和球坐标下的 ...
Matlab可视化四维数据（用颜色表示柱坐标上某点的数值大小）
Matlab新手,只是粗略看过一遍b站的教程.在使用Matlab的过程中,深感函数使用.编程思维等的困难.感觉还是要加强编程思维的训练比较好. 书归正传,这是我第一次运用Matlab绘制三维图.需要将 ...
柱坐标系下的ns方程_麦克斯韦方程组小结
一.▽ 算子.点积.叉积 l▽ 算子叫"del"算子,即<< span="">∂/∂x,∂/∂y,∂/∂z>,可以理解为一个符号向量,向 ...
多重网格法-松弛迭代法-二维泊松方程-python实现
这几天在家躲避疫情,闲来无事,写了这个多重网格法求解泊松方程的算法的代码. 多重网格法可能是目前为止解泊松方程最快的算法,n个格点需要n次计算就可以收敛,而快速傅里叶变换的收敛速度是n*logn, 共 ...
多重网格法（multigrid）求解1d泊松方程--python
多重网格法(multigrid)求解1d泊松方程多重网格法可能是目前已知最快的求解泊松方程的算法,特别是在边界复杂的情况下,实用性远高于快速傅里叶算法. 我写里写了一个一维求解泊松方程的程序,供大家 ...
AidLux“换脸”案例源码详解（Python）
"换脸"案例源码详解 (Python) faceswap_gui.py用于换脸,可与facemovie_gui.py身体互换源码(上一篇文章)对照观看打开faceswap_gui ...
利用python处理dna序列_详解基于python的全局与局部序列比对的实现(DNA)
程序能实现什么 a.完成gap值的自定义输入以及两条需比对序列的输入 b.完成得分矩阵的计算及输出 c.输出序列比对结果 d.使用matplotlib对得分矩阵路径的绘制一.实现步骤 1.用户输入步 ...

柱坐标下多重网格法解泊松方程-python

柱坐标下多重网格法解泊松方程-python相关推荐

最新文章

热门文章