burgers equation 1.1公式

burgers equation 1.1公式

2D

公式1(python):

laplace_operator

def apply_laplacian(mat, dx = 1.0):# dx is inversely proportional to N"""This function applies a discretized Laplacianin periodic boundary conditions to a matrixFor more information see """# the cell appears 4 times in the formula to compute# the total differenceneigh_mat = -py()# Each direct neighbor on the lattice is counted in# the discrete difference formulaneighbors = [ ( 4/3,  (-1, 0) ),( 4/3,  ( 0,-1) ),( 4/3,  ( 0, 1) ),( 4/3,  ( 1, 0) ),(-1/12,  (-2, 0)),(-1/12,  (0, -2)),(-1/12,  (0, 2)),(-1/12,  (2, 0)),]# shift matrix according to demanded neighbors# and add to this cell with corresponding weightfor weight, neigh in neighbors:neigh_mat += weight * np.roll(mat, neigh, (0,1))return neigh_mat/dx**2

ll()

np.roll()

mat=[0, 1, 2, 3, 4,5, 6, 7, 8, 9,4,7]
maat&#shape(mat, (3,4))
print(maat)
neigh=(-1,0)
a&#ll(maat, neigh, (1,0))
print(a)

result

[[0 1 2 3][4 5 6 7][8 9 4 7]]
[[1 2 3 0][5 6 7 4][9 4 7 8]]

分析因为 n e i g h = ( − 1 , 0 ) neigh=(-1,0) neigh=(−1,0),向左一格

apply_dx

def apply_dx(mat, dx = 1.0):''' central diff for dx'''# np.roll, axis=0 -> row # the total differenceneigh_mat = -py()# Each direct neighbor on the lattice is counted in# the discrete difference formulaneighbors = [ ( 1.0/12,  (2, 0) ),( -8.0/12,  (1, 0) ),( 8.0/12,  (-1, 0) ),( -1.0/12,  (-2, 0) )]# shift matrix according to demanded neighbors# and add to this cell with corresponding weightfor weight, neigh in neighbors:neigh_mat += weight * np.roll(mat, neigh, (0,1))return neigh_mat/dx

apply_dy

def apply_dy(mat, dy = 1.0):''' central diff for dx'''# the total differenceneigh_mat = -py()# Each direct neighbor on the lattice is counted in# the discrete difference formulaneighbors = [ ( 1.0/12,  (0, 2) ),( -8.0/12,  (0, 1) ),( 8.0/12,  (0, -1) ),( -1.0/12,  (0, -2) )]# shift matrix according to demanded neighbors# and add to this cell with corresponding weightfor weight, neigh in neighbors:neigh_mat += weight * np.roll(mat, neigh, (0,1))return neigh_mat/dy

get_temporal_diff

def get_temporal_diff(U, V, R, dx):# u and v in (h, w)laplace_u = apply_laplacian(U, dx)laplace_v = apply_laplacian(V, dx)u_x = apply_dx(U, dx)v_x = apply_dx(V, dx)u_y = apply_dy(U, dx)v_y = apply_dy(V, dx)# governing equationu_t = (1.0/R) * laplace_u - U * u_x - V * u_yv_t = (1.0/R) * laplace_v - U * v_x - V * v_yreturn u_t, v_t

update_rk4

def update_rk4(U0, V0, R=100, dt=0.05, dx=1.0):"""Update with Runge-kutta-4 methodSee """############# Stage 1 ############### compute the diffusion part of the updateu_t, v_t = get_temporal_diff(U0, V0, R, dx)K1_u = u_tK1_v = v_t############# Stage 2 ##############U1 = U0 + K1_u * dt/2.0V1 = V0 + K1_v * dt/2.0u_t, v_t = get_temporal_diff(U1, V1, R, dx)K2_u = u_tK2_v = v_t############# Stage 3 ##############U2 = U0 + K2_u * dt/2.0V2 = V0 + K2_v * dt/2.0u_t, v_t = get_temporal_diff(U2, V2, R, dx)K3_u = u_tK3_v = v_t############# Stage 4 ##############U3 = U0 + K3_u * dtV3 = V0 + K3_v * dtu_t, v_t = get_temporal_diff(U3, V3, R, dx)K4_u = u_tK4_v = v_t# Final solutionU = U0 + dt*(K1_u+2*K2_u+2*K3_u+K4_u)/6.0V = V0 + dt*(K1_v+2*K2_v+2*K3_v+K4_v)/6.0return U, V

main

part1

    M, N = 256, 256n_simu_steps = 30000dt = 0.003 # maximum 0.003dx = 1.0 / MR = 200.0# get initial condition from random fielddevice = torch.device('cuda')GRF = GaussianRF(2, M, alpha=2, tau=5, device=device)U, V = GRF.sample(2) # U and V have shape of [128, 128]U = U.cpu().numpy()V = V.cpu().numpy()U_record = U.copy()[None,...]V_record = V.copy()[None,...]

note:

1. 此处利用random_field得到U,V初始值

U, V = GRF.sample(2)

得到的GRF.sample(2)返回值 u u u
s i z e = [ 2 , 256 , 256 ] size=[2,256,256] size=[2,256,256]
U = u [ 0 , : , : ] U=u[0,:,:] U=u[0,:,:]
V = u [ 1 , : , : ] V=u[1,:,:] V=u[1,:,:]

2. U_record、V_record维数增加
   利用[None,…]

U_record = U.copy()[None,...]
V_record = V.copy()[None,...]

part2
循环计算

    for step in range(n_simu_steps):U, V = update_rk4(U, V, R, dt, dx) #[h, w]if (step+1) % 20 == 0:print(step, 'n')U_record = np.concatenate((U_record, U[None,...]), axis=0) # [t,h,w]V_record = np.concatenate((V_record, V[None,...]), axis=0)UV = np.concatenate((U_record[None,...], V_record[None,...]), axis=0) # (c,t,h,w)UV = np.transpose(UV, [1, 0, 2, 3]) # (t,c,h,w)

note:

默认按照行数叠加矩阵

np.concatenate()

交换维数

np.transpose()

random_field

与原版本不同，因为得用cuda版本，而我cuda版本太低所以换了计算方法
原版本
↓ downarrow ↓

import torch
import math
from torch.fft import ifftimport matplotlib.pyplot as plt
import matplotlibfrom timeit import default_timer
#device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
class GaussianRF(object):def __init__(self, dim, size, alpha=2, tau=3, sigma=None, boundary="periodic", device=None):#device = torch.device("cuda" if torch.cuda.is_available() else "cpu")self.dim = dimself.device = deviceif sigma is None:sigma = tau**(0.5*(2*alpha - self.dim))k_max = size//2if dim == 1:k = torch.cat((torch.arange(start=0, end=k_max, step=1, device=device), torch.arange(start=-k_max, end=0, step=1, device=device)), 0)self.sqrt_eig = size*math.sqrt(2.0)*sigma*((4*(math.pi**2)*(k**2) + tau**2)**(-alpha/2.0))self.sqrt_eig[0] = 0.0elif dim == 2:wavenumers = torch.cat((torch.arange(start=0, end=k_max, step=1, device=device), torch.arange(start=-k_max, end=0, step=1, device=device)), 0).repeat(size,1)k_x = anspose(0,1)k_y = wavenumersself.sqrt_eig = (size**2)*math.sqrt(2.0)*sigma*((4*(math.pi**2)*(k_x**2 + k_y**2) + tau**2)**(-alpha/2.0))self.sqrt_eig[0,0] = 0.0elif dim == 3:wavenumers = torch.cat((torch.arange(start=0, end=k_max, step=1, device=device), torch.arange(start=-k_max, end=0, step=1, device=device)), 0).repeat(size,size,1)k_x = anspose(1,2)k_y = wavenumersk_z = anspose(0,2)self.sqrt_eig = (size**3)*math.sqrt(2.0)*sigma*((4*(math.pi**2)*(k_x**2 + k_y**2 + k_z**2) + tau**2)**(-alpha/2.0))self.sqrt_eig[0,0,0] = 0.0self.size = []for j in range(self.dim):self.size.append(size)self.size = tuple(self.size)def sample(self, N):coeff = torch.randn(N, *self.size, 2, device=self.device)coeff[...,0] = self.sqrt_eig*coeff[...,0]coeff[...,1] = self.sqrt_eig*coeff[...,1]# u = torch.fft.ifft(coeff, self.dim, normalized=False)u = torch.fft.ifft(coeff, self.dim)u = u[...,0]return u

改变后
↓ downarrow ↓

import torch
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlibfrom timeit import default_timer
import scipy.io as scioclass GaussianRF(object):def __init__(self, dim, size, alpha=2, tau=3, sigma=None, boundary="periodic", device=None):self.dim = dimself.device = deviceif sigma is None:sigma = tau**(0.5*(2*alpha - self.dim))k_max = size//2if dim == 1:k = torch.cat((torch.arange(start=0, end=k_max, step=1, device=device), torch.arange(start=-k_max, end=0, step=1, device=device)), 0)self.sqrt_eig = size*math.sqrt(2.0)*sigma*((4*(math.pi**2)*(k**2) + tau**2)**(-alpha/2.0))self.sqrt_eig[0] = 0.0elif dim == 2:wavenumers = torch.cat((torch.arange(start=0, end=k_max, step=1, device=device), torch.arange(start=-k_max, end=0, step=1, device=device)), 0).repeat(size,1)k_x = anspose(0,1)k_y = wavenumersdata_dir = 'C:/Users/admin/Documents/MATLAB/Examples/pde/pdedemo5/co.mat'data = scio.loadmat(data_dir)self.sqrt_eig = data['sqrt_eig']self.sqrt_eig[0,0] = 0.0elif dim == 3:wavenumers = torch.cat((torch.arange(start=0, end=k_max, step=1, device=device), torch.arange(start=-k_max, end=0, step=1, device=device)), 0).repeat(size,size,1)k_x = anspose(1,2)k_y = wavenumersk_z = anspose(0,2)self.sqrt_eig = (size**3)*math.sqrt(2.0)*sigma*((4*(math.pi**2)*(k_x**2 + k_y**2 + k_z**2) + tau**2)**(-alpha/2.0))self.sqrt_eig[0,0,0] = 0.0self.size = []for j in range(self.dim):self.size.append(size)self.size = tuple(self.size)def sample(self, N):# coeff = torch.randn(N, *self.size, 2, device=self.device)## coeff[...,0] = self.sqrt_eig*coeff[...,0]# coeff[...,1] = self.sqrt_eig*coeff[...,1]coeff = torch.from_numpy(np.random.randint( size=(2, 256, 256, 2))).double()# coeff&#(device)# [t,c,h,w]print(self.sqrt_eig)print(self.sqrt_eig.shape)print(coeff[:, :, :, 0] .shape)coeff[:, :, :, 0] = torch.from_numpy(self.sqrt_eig) * coeff[:, :, :, 0]coeff[..., 1] = torch.from_numpy(self.sqrt_eig) * coeff[:, :, :, 1]u = torch.ifft(coeff, self.dim, normalized=False)u = u[...,0]return u

加入一些print()方便观察矩阵维度
加的是利用matlab计算的后续进行matlab的分析