Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- nx=3
- nz=4
- flux2 = np.empty((nx+1, nz+1))
- flux2[:] = np.random.normal(size=(nx+1, nz+1))
- flux3 = flux2.copy()
- bz = np.random.normal(size=(nx+1, nz+1))
- bx = np.random.normal(size=(nx+1, nz+1))
- dz = 0.112344
- dx = 0.015543
- for i in np.arange(1,nx+1):
- for j in np.arange(1,nz+1):
- flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)
- a = 0.5
- aexp = np.arange(nz).reshape(nz, 1) - np.arange(nz).reshape(1, nz)
- abcoeff = 0.5**aexp
- abcoeff[aexp<0] = 0
- for i in np.arange(1,nx+1):
- b = 0.5*flux3[i-1, 1:] + 0.5*bz[i-1, 1:]*dx - 0.5*bx[i,:-1]*dz
- bvals = (abcoeff * b.reshape(1, nz)).sum(axis=1)
- n = np.arange(1, nz+1)
- x0 = flux3[i, 0]
- #ab = b * np.cumsum(a**(n-1))
- flux3[i, 1:] = a**n * x0 + bvals
- print np.abs(flux2-flux3) < 1e-15
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement