Commit 0d29205b authored by Nikolai G. Lehtinen's avatar Nikolai G. Lehtinen
Browse files

initial commit

parents
**~
__pycache__/
trump/__pycache__/
TRUMP (Total Removal of Unwanted Matrix Programming)
=======
# Introduction
TRUMP (Total Removal of Unwanted Matrix Programming) is a package whose goal is
to help people with programming various finite-difference schemes. It allows the
user
1. to store the discretized values for physical variables in an object called
_Field_ which is defined on a 1D or 2D equally-spaced rectangular grid with
ghost cells at boundaries which are accessible with negative indices or indices
which look like `end+1` etc. The intermediate results for operations with fields
are stored in objects of class _ExtendedArray_ which have the same indexing features,
(but cannot generate operators or boundary conditions, see the rest of items).
2. to create finite-difference operators which act on Fields;
3. to specify boundary conditions on the Field which will be either
* easily updated when using forward finite-differencing schemes
* automatically enforced when inverting finite-difference operators acting on this field;
4. to set up a solver which will invert the finite-difference operator.
The purpose of the package is obviously to make finite-difference coding great again.
# Getting started: Running examples
Please use **Python3** (any version). Other requirements are listed in the [`requirements.txt`](requirements.txt) file.
We emphasize that the package does not provide any numerical solvers, it only acts as an interface
to solvers in Numpy. It is essential to have a good version of NumPy, because some of the previous versions were buggy.
## Read the manual
Read the [manual](trump.pdf) (at least, the introduction) and run [`trump_intro.py`](trump_intro.py)
## Operators, boundary conditions and solvers
You can start with running [`trump_demo.py`](trump_demo.py). The file is very short, please read what it does.
## Advection equation solution
Then, more complicated examples which demonstrate solution of advection equation
```math
\frac{\partial{f}}{\partial{t}} = -\nabla\cdot(\mathbf{v}f)
```
We use several different methods to solve it:
* First order - CIR method
* Second order - MacCormack method (a more stable modification of Lax-Wendroff)
* Third order - Leonard method
The featured third-order method is of upwind type and also uses a flux corrector due to the same
author. To quote the author, "the fourth-order diffusion kills the third-order dispersion", therefore
the results are even better than if the fourth-order method was used.
The files are:
* [`car_advection_test.py`](car_advection_test.py) - rotation in cartesian coordinates. It uses
divergenceless velocity given by
```math
\mathbf{v} = \omega(\hat{\mathbf{z}}\times\mathbf{r})
```
* [`cyl_advection_test.py`](cyl_advection_test.py) - motion in cylindrical coordinates (there is an
option for divergenceless velocity, too)
* A more advanced example: [`cyl_advection_symmetry_test.py`](more_examples/cyl_advection_symmetry_test.py) - Here we test
a situation for which the density changes uniformly, to make sure that
behavior at the axis is correct.
The former location of the project was [here](https://gitlab.uib.no/BCSS-Q4/TRUMP).
# Other examples of _Markdown_ syntax
Gitlab uses [*GitLab Flavored Markdown (GFM)*](https://gitlab.com/help/user/markdown).
(There is another [copy of that manual](https://docs.gitlab.com/ee/user/markdown.html) which
is not rendered very well). [Cheatsheet](https://guides.github.com/pdfs/markdown-cheatsheet-online.pdf).
[For developers](https://github.github.com/gfm/) (don't read if you don't need problems).
Original flavor (?) [MarkDown syntax](https://daringfireball.net/projects/markdown/syntax) and
[cheatsheet](https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet).
Images: ![Image](velocity.png)
Block comments: use `<!--` and `-->`:
<!--
All this should be
commented out
-->
<!-- and this too -->
Another way: `[//] (comment)` (parentheses are **obligatory**!).
[//]: # (should be invisible)
Inline comments: use a trick of an empty link: `[](comment here)` [](comment here).
Math $`a^2+b^2=c^2`$
Two spaces at the end of a line
leave a line break.
Horizontal rule:
---
> Markdown uses email-style > characters for blockquoting.
Inline <abbr title="Hypertext Markup Language">HTML</abbr> is supported.
This was ~~useless~~ useful [:smiley:](https://www.webpagefx.com/tools/emoji-cheat-sheet/)!
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Dec 4 19:18:10 2016
Solution of advective equation:
df/dt + div(v*f) = 0
using various linear schemes.
f in div(v*f) is represented as a "face value" of f
Literature:
Leonard [1979] - 3rd order upwind scheme in 1D (called QUICKEST)
Leonard and Niknafs [1991] - in 2D (called UTOPIA)
Leonard [1991] - the flux-limiting algorithm in 1D (called ULTIMATE)
Previous versions:
advection2d_rotation_symbolic.py
advection2d_rotation_waverchoice.py
advection2d_rotation_demo.py
@author: nle003
"""
#%% Preliminaries
# Imports
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from my_utils import myfloat
import time
import trump
trump.DEBUGLEVEL=0
trump.SET_SPARSE(True)
from trump import ExtendedArray2D, end, Grid, grid_x,\
c_dif, n_dif, c_ave, n_ave, c_upc, n_upc, Field2D
from trump_2D import Geometry2D_car, AdvectionStepI, UpdateablePcolor
from trump_2D import c_difx, c_dify, n_difx, n_dify
def n_div(fx,fy): return n_dif(fx,0)+n_dif(fy,1)
def c_upx(f,v): return c_upc(f,v,0)
def c_upy(f,v): return c_upc(f,v,1)
################################################################################
#%% Testing
vmax = 0.45
do_opposite_rotation=False
num_smear = 0
do_log_plot = False
# Setup the grid and velocity field: CCW for r<R1, CW for R1<r<R2, w=vmax/R2
# The unknown function f is calculated on a Nx x Ny array, at x,y coordinates
Lx = 2; Ly = 2; dx=.01; dy=.01;
#Lx = 200; Ly=200; dx=1; dy=1
Nx=np.int(Lx/dx)
Ny=np.int(Ly/dy)
gridx = Grid(n_cells=Nx,delta=dx,start=-Lx/2)
gridy = Grid(n_cells=Ny,delta=dy,start=-Ly/2)
# Extended coordinates are for pcolor plots
#xe = geom_x(geomx,0)
#ye = geom_x(geomy,0)
#x=geom_x(geomx,1)
#y=geom_x(geomy,1)
# Velocity (scaled), calculated on the grid edges as v = curl (psi z)
#R2=np.min([Lx/2-3*dx,Ly/2-3*dy])
R2 = np.sqrt(2)*Lx/2
if do_opposite_rotation:
R1 = np.min([Lx/4,Ly/4])
else:
R1 = 0
vangle=vmax/R2 # angular velocity
# Zero-divergence velocity, curl of psi*z_unit
geom = Geometry2D_car(gridx,gridy,nls=(3,3),nus=(3,3))
psi=ExtendedArray2D(geom,stags=(1,1),nls=(3,3),nus=(3,3))
yg,xg=np.meshgrid(psi.xe(1),psi.xe(0))
rg=np.sqrt(xg**2+yg**2)
psi.arr = -rg**2/2*vangle # at points (i-1/2,j-1/2)
i2=(rg>=R1)
psi.arr[i2]=rg[i2]**2/2*vangle-R1**2*vangle
i3=(rg>=R2)
psi.arr[i3]=R2**2/2*vangle-R1**2*vangle
vx = n_dify(psi) # vx at points (i-1/2,j)
vy = -n_difx(psi) # vy at points (i,j-1/2)
# Check if the divergence is zero at points (i,j)
print('max div v =',np.max(np.abs(n_div(vx,vy)).arr),flush=True)
dt = 1/(np.max(vx.arr)/dx+np.max(vy.arr)/dy)
################################################################################
#%% Prepare the figures
methods=['CIR','MacCormack','MacCormack_Reverse','UTOPIA','ULTIMATE']
nfig=len(methods)
savedirs=[None]*nfig
info = lambda n,f,m: '[{:.3f},{:.3f}],n={:.2f},{:s}'.format(np.min(f.arr),np.max(f.arr),n,m)
zmin,zmax = (-8,.5) if do_log_plot else (-.4,1.4)
fig=[UpdateablePcolor(100+kfig,figsize=(6,5),savedir=savedirs[kfig],cmap=mpl.cm.jet,zmin=zmin,zmax=zmax)
for kfig in range(nfig)]
#%%#############################################################################
# Initial value
fini = ExtendedArray2D(geom,stags=(0,0),nls=(2,2),nus=(2,2))
fini.alloc(np.double)
yg,xg=np.meshgrid(fini.xe(1),fini.xe(0))
xsqc=0.; ysqc=Ly/6.
ax=Lx/4; ay=Ly/6
fini.arr[(np.abs(xg-xsqc)<ax) & (np.abs(yg-ysqc)<ay)]=1.
def smear(ea,n=0):
for k in range(n):
ea.setv = n_ave(c_ave(n_ave(c_ave(ea,0),0),1),1)
#return ea
smear(fini,num_smear)
################################################################################
#%% The main cycle
# Initial value
farr=[]
for kfig in range(nfig):
if kfig>=1: dd = 2
else: dd = 1
f0=Field2D('f'+str(kfig+1),geom,stags=(0,0),nls=(dd,dd),nus=(dd,dd))
tt1 = time.time()
if dd==1:
f0.bc[-1,-1:end+1]=0
f0.bc[end+1,-1:end+1]=0
f0.bc[:,-1]=0
f0.bc[:,end+1]=f0.bc[:,end]
else:
f0.bc[-2:-1,-2:end+2]=0
f0.bc[end+1:end+2,-2:end+2]=0
f0.bc[:,-2:-1]=0
f0.bc[:,end+1]=f0.bc[:,end]
f0.bc[:,end+2]=f0.bc[:,end]
tt2 = time.time()
f0.bc.freeze() # takes pretty long
tt3 = time.time()
print('To record: t = ',tt2-tt1,', to freeze = ',tt3-tt2,flush=True)
farr.append(f0)
#ye = f0.pcolor_x(1)
for f0 in farr:
f0.setv = fini
totals = [geom.integrate(f) for f in farr]
print('totals =',totals)#xe = f0.pcolor_x(0)
#%%
for kt in range(50001):
for kfig in range(nfig):
farr[kfig] += AdvectionStepI(farr[kfig],vx*dt,vy*dt,methods[kfig])
farr[kfig].bc.apply()
#f = farr[-1]
#fv1, dfv3 = FaceValue_UTOPIA(f,vx*dt,vy*dt)
#fluxx,fluxy = TVD_ULTIMATE(f,(vx*dt,vy*dt),fv1,dfv3,apply=True)
#f -= geom.Div(fluxx,fluxy)
#f.bc.apply()
if kt<100 or (kt<500 and kt % 10==0) or kt % 100==0:
n_turns = vangle*dt*kt/(2*np.pi)
for kfig in range(nfig):
f = farr[kfig]
dtot = geom.integrate(f)-totals[kfig]
fplot = np.log10(f+1e-8) if do_log_plot else f
fig[kfig].plot(fplot,title=info(n_turns,f,methods[kfig]+' '+str(myfloat(dtot))))
plt.pause(0.1)
# fig[0].wait() # progress by clicking mouse in Figure 100
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Nov 18 18:48:23 2017
Test of advective schemes in cylindrical coordinates.
The velocity profile we test is (such that div v = 0, vz=0 at z=+-L/2 and vr=0 at r=R):
vr = vo (π/L) J_1(br/R) sin(πz/L) = -dψ/dz
vz = vo bJ_0(br/R) cos(πz/L) = (1/r)d(rψ)/dr
with b being the first root of J_1, b=spec.jn_zeros(1,1)[0]=3.8317
It may be represented as curl A where A_φ = ψ and is given by
ψ(r,z) = vo J_1(br/R) cos(πz/L)
@author: nle003
"""
#%% Preliminaries
import time
from my_utils import np,plt,spec,mpl,myfloat
π = np.pi
b=spec.jn_zeros(1,1)[0]
import trump
trump.DEBUGLEVEL=2
from trump import end, span, Grid, ExtendedArray2D, Field2D, c_ave, n_dif
from trump_2D import Geometry2D_cyl, UpdateablePcolor, AdvectionStepI
#def c_avez(f): return c_ave(f,1)
#%% Geometry setup
vo = 1
R = 1
L = 2
Nr = 100
Nz = 200
dr = R/Nr
dz = L/Nz
gridr = Grid(n_cells=Nr,delta=dr,start=0)
gridz = Grid(n_cells=Nz,delta=dz,start=-L/2)
cyl = Geometry2D_cyl(gridr,gridz,nls=(2,2),nus=(2,2))
#%% Velocity
divergence_free=True
if divergence_free:
# Zero-divergence velocity, curl of psi*phi_unit
psi=ExtendedArray2D(cyl,stags=(1,1),nls=(3,3),nus=(3,3))
zpg,rpg=np.meshgrid(psi.xe(1),psi.xe(0))
if False:
b=spec.jn_zeros(1,1)[0]
psi.arr=vo*spec.j1(b*rpg/R)*np.cos(π*zpg/L)
else:
# Psi can actually be arbitrary, just make sure vr=0 at r=0
psi.arr = vo*np.sin(π*rpg/R)*np.cos(π*zpg/L)
vr = -n_dif(psi,1) # vx at points (i-1/2,j)
vz = n_dif(psi*c_ave(cyl.rd,1),0)/c_ave(cyl.ri,1) # vy at points (i,j-1/2)
else:
nls_r=cyl.r.nls
vr = ExtendedArray2D(cyl,stags=(1,0),nls=cyl.rd.nls,nus=cyl.ri.nus)
zvg,rvg = np.meshgrid(vr.xe(1),vr.xe(0))
vr.arr = vo*np.sin(π*rvg/R)*np.sin(π*zvg/L)
vz = ExtendedArray2D(cyl,stags=(0,1),nls=cyl.ri.nls,nus=cyl.ri.nus)
zvg,rvg = np.meshgrid(vz.xe(1),vz.xe(0))
vz.arr = vo*np.cos(π*rvg/R)*np.cos(π*zvg/L)
# Check if the divergence is zero at points (i,j)
print('max div v =',np.max(np.abs(cyl.Div(vr,vz)).arr),flush=True)
dt = min(dr,dz)/max(np.max(vr.arr),np.max(vz.arr))
#%% Initial value
methods = ['CIR','MacCormack','MacCormack_Reverse','UTOPIA','ULTIMATE']
nfig = len(methods)
fs=[]
totals = np.zeros((nfig,))
if False:
rtmp = cyl.ri.copy()
rtmp.arr=cyl.ri.arr.copy()
rtmp[0,span]=dr/8
def get_tot(f):
return 2*np.pi*dz*dr*np.sum((f*rtmp)[:,:])
else:
def get_tot(f):
return cyl.integrate(f)
for kfig in range(nfig):
if methods[kfig] in ['CIR']:
dd=1
elif methods[kfig] in ['MacCormack','MacCormack_Reverse','UTOPIA','ULTIMATE']:
dd=2
else:
raise Exception('Unknown method')
f0 = Field2D('f'+str(kfig+1),cyl,stags=(0,0),nls=(dd+1,dd),nus=(dd,dd))
tt1 = time.time()
for k in range(dd+1):
f0.bc[-k-1,:end-1]=f0.bc[k+1,:end-1] # The values at nz=end are dependent due to periodic bc
f0.bc[:,end+k]=f0.bc[:,k]
for k in range(1,dd+1):
f0.bc[:,-k]=f0.bc[:,end-k]
tt2 = time.time()
assert f0.bc.num>0 # we have set at least some boundary conditions
f0.bc.freeze() # takes pretty long
tt3 = time.time()
print('To record: t = ',tt2-tt1,', to freeze = ',tt3-tt2,flush=True)
zfg,rfg=np.meshgrid(f0.xe(1),f0.xe(0))
zc=0;
ar=2*R/3; az=L/4
f0.arr[(np.abs(zfg-zc)<az) & (np.abs(rfg)<ar)]=1.
totals[kfig]=get_tot(f0)
fs.append(f0)
print('totals =',totals)
#%%
fig = [UpdateablePcolor(100+k,figsize=(6,5),zmin=-0.4,zmax=1.4,cmap=mpl.cm.jet)
for k in range(nfig)]
info = lambda k,f,m: '[{:.3g},{:.3g}],k={},{:s}'.format(np.min(f.arr),np.max(f.arr),k,m)
rc,zc = cyl.xcs
#%%
for kt in range(201):
for kfig in range(nfig):
f = fs[kfig]
#tmp = cyl.AdvectionStepI(f,vr*dt,vz*dt,methods[kfig])
#print(get_tot(tmp))
f += AdvectionStepI(f,vr*dt,vz*dt,methods[kfig])
#f[f<0]=0
f.bc.apply()
#dtot = cyl.integrate(f,rb=np.array([0,rc[end]]),zb=np.array([zc[0],zc[end]]))-totals[kfig]
dtot = get_tot(f)-totals[kfig]
if kt<100 or (kt<500 and kt % 10==0) or kt % 100==0:
for kfig in range(nfig):
fig[kfig].plot(fs[kfig],title=info(kt,fs[kfig],methods[kfig]+' '+str(myfloat(dtot))))
plt.pause(0.1)
# fig[0].wait() # progress by clicking mouse in Figure 100
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Nov 18 18:48:23 2017
Test of advective schemes in cylindrical coordinates.
This particular file is to dispel any doubts that velocity on axis is handled correctly. And it is!
The velocity profile we test is (such that div v = 0, vz=0 at z=+-L/2 and vr=0 at r=R):
vr = vo (π/L) J_1(br/R) sin(πz/L) = -dψ/dz
vz = vo bJ_0(br/R) cos(πz/L) = (1/r)d(rψ)/dr
with b being the first root of J_1, b=spec.jn_zeros(1,1)[0]=3.8317
It may be represented as curl A where A_φ = ψ and is given by
ψ(r,z) = vo J_1(br/R) cos(πz/L)
@author: nle003
"""
#%% Preliminaries
import time
from my_utils import np,plt,spec,mpl,myfloat
π = np.pi
b=spec.jn_zeros(1,1)[0]
import trump
trump.DEBUGLEVEL=2
from trump import end, span, Grid, ExtendedArray2D, Field2D, c_ave, n_dif
from trump_2D import Geometry2D_cyl, UpdateablePcolor, AdvectionStepI
#def c_avez(f): return c_ave(f,1)
#%% Geometry setup
vo = 1
R = 1
L = 2
Nr = 50
Nz = 100
dr = R/Nr
dz = L/Nz
gridr = Grid(n_cells=Nr,delta=dr,start=0)
gridz = Grid(n_cells=Nz,delta=dz,start=-L/2)
cyl = Geometry2D_cyl(gridr,gridz,nls=(2,2),nus=(2,2))
#%% Velocity
def vabs(r):
return r
vr = ExtendedArray2D(cyl,stags=(1,0),nls=cyl.rd.nls,nus=cyl.ri.nus)
zvg,rvg = np.meshgrid(vr.xe(1),vr.xe(0))
rabs = np.abs(rvg**2+zvg**2)
vr.arr = vabs(rabs)*rvg/rabs
vz = ExtendedArray2D(cyl,stags=(0,1),nls=cyl.ri.nls,nus=cyl.ri.nus)
zvg,rvg = np.meshgrid(vz.xe(1),vz.xe(0))
rabs = np.abs(rvg**2+zvg**2)
vz.arr = -2*vabs(rabs)*zvg/rabs
# Check if the divergence is zero at points (i,j)
print('max div v =',np.max(np.abs(cyl.Div(vr,vz)).arr),flush=True)
dt = min(dr,dz)/(np.max(np.abs(vr.arr))+np.max(np.abs(vz.arr)))
#%% Initial value
methods = ['CIR','MacCormack','MacCormack_Reverse','UTOPIA','ULTIMATE']
nfig = len(methods)
fs=[]
totals = np.zeros((nfig,))
if False:
rtmp = cyl.ri.copy()
rtmp.arr=cyl.ri.arr.copy()
rtmp[0,span]=dr/8
def get_tot(f):
return 2*np.pi*dz*dr*np.sum((f*rtmp)[:,:])
else:
def get_tot(f):
return cyl.integrate(f)
for kfig in range(nfig):
if methods[kfig] in ['CIR']:
dd=1
elif methods[kfig] in ['MacCormack','MacCormack_Reverse','UTOPIA','ULTIMATE']:
dd=2
else:
raise Exception('Unknown method')
f0 = Field2D('f'+str(kfig+1),cyl,stags=(0,0),nls=(dd+1,dd),nus=(dd,dd))
tt1 = time.time()
for k in range(dd+1):
f0.bc[-k-1,:end-1]=f0.bc[k+1,:end-1] # The values at nz=end are dependent due to periodic bc
f0.bc[:,end+k]=f0.bc[:,k]
for k in range(1,dd+1):
f0.bc[:,-k]=f0.bc[:,end-k]
tt2 = time.time()
assert f0.bc.num>0 # we have set at least some boundary conditions
f0.bc.freeze() # takes pretty long
tt3 = time.time()
print('To record: t = ',tt2-tt1,', to freeze = ',tt3-tt2,flush=True)
zfg,rfg=np.meshgrid(f0.xe(1),f0.xe(0))
zc=0;
ar=R/2; az=L/4
f0.arr[(np.abs(zfg-zc)<az) & (np.abs(rfg)<ar)]=1.
totals[kfig]=get_tot(f0)
fs.append(f0)
print('totals =',totals)
#%%
fig = [UpdateablePcolor(100+k,figsize=(6,5),zmin=-0.4,zmax=1.4,cmap=mpl.cm.jet)
for k in range(nfig)]
info = lambda k,f,m: '[{:.3g},{:.3g}],k={},{:s}'.format(np.min(f.arr),np.max(f.arr),k,m)
rc,zc = cyl.xcs
#%%
for kt in range(201):
for kfig in range(nfig):
f = fs[kfig]
#tmp = cyl.AdvectionStepI(f,vr*dt,vz*dt,methods[kfig])
#print(get_tot(tmp))
f += AdvectionStepI(f,vr*dt,vz*dt,methods[kfig])
#f[f<0]=0
f.bc.apply()
#dtot = cyl.integrate(f,rb=np.array([0,rc[end]]),zb=np.array([zc[0],zc[end]]))-totals[kfig]
dtot = get_tot(f)-totals[kfig]
if kt<100 or (kt<500 and kt % 10==0) or kt % 100==0:
for kfig in range(nfig):
fig[kfig].plot(fs[kfig],title=info(kt,fs[kfig],methods[kfig]+' '+str(myfloat(dtot))))
plt.pause(0.1)
# fig[0].wait() # progress by clicking mouse in Figure 100
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 17 10:33:45 2017
Some utilities of questionable usability (some may be reinventions of a wheel)
@author: nle003
"""
import pickle
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import io
import scipy.constants as phys
import scipy.special as spec
import scipy.interpolate
####################################################################################################
#%% A few general-purpose routines
def line_interp2d(x,y,xp,yp,zp):
"""Workaround around clumsy 2D interpolation in scipy. Stolen from
https://stackoverflow.com/questions/47087109/
evaluate-the-output-from-scipy-2d-interpolation-along-a-curve"""
shape=x.shape
assert y.shape==shape
f = scipy.interpolate.interp2d(xp, yp, zp, bounds_error=True)
return scipy.interpolate.dfitpack.bispeu(
f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4],
x.flat, y.flat)[0].reshape(shape)
def ones_as(x):
"Scalar or vector ones"
return np.ones(np.asarray(x).shape)[()]
def zeros_as(x):
"Scalar or vector zeros"
return np.zeros(np.asarray(x).shape)[()]