Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

work in progress #6

Open
wants to merge 8 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -23,4 +23,11 @@ docs/site/
# environment.
Manifest.toml

.vscode
.vscode

#Julia system image
JuliaSysimage.dll

#python stuff
__pycache__
*.pyc
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ version = "0.1.2"

[deps]
ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -15,3 +16,4 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
Binary file not shown.
352 changes: 352 additions & 0 deletions Python Burgers example/benchmarkedBurgers.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,352 @@
# 2nd-order accurate finite-volume implementation of the inviscid Burger's
# equation with piecewise linear slope reconstruction
#
# We are solving u_t + u u_x = 0 with outflow boundary conditions
#
# M. Zingale (2013-03-26)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import sys
import time
import cProfile
from memory_profiler import profile
from pympler import asizeof
from pympler import muppy
from pympler import summary
from pympler import classtracker
import tracemalloc
import linecache
import os

mpl.rcParams['mathtext.fontset'] = 'cm'
mpl.rcParams['mathtext.rm'] = 'serif'

mpl.rcParams['font.size'] = 12
mpl.rcParams['legend.fontsize'] = 'large'
mpl.rcParams['figure.titlesize'] = 'medium'

class Grid1d(object):
def __init__(self, nx, ng, xmin=0.0, xmax=1.0, bc="outflow"):

self.nx = nx
self.ng = ng

self.xmin = xmin
self.xmax = xmax

self.bc=bc

# python is zero-based. Make easy intergers to know where the
# real data lives
self.ilo = ng
self.ihi = ng+nx-1

# physical coords -- cell-centered, left and right edges
self.dx = (xmax - xmin)/(nx)
self.x = xmin + (np.arange(nx+2*ng)-ng+0.5)*self.dx

# storage for the solution
self.u = np.zeros((nx+2*ng), dtype=np.float64)

def scratch_array(self):
""" return a scratch array dimensioned for our grid """
return np.zeros((self.nx+2*self.ng), dtype=np.float64)

def fill_BCs(self):
""" fill all ghostcells as periodic """

if self.bc == "periodic":

# left boundary
self.u[0:self.ilo] = self.u[self.ihi-self.ng+1:self.ihi+1]

# right boundary
self.u[self.ihi+1:] = self.u[self.ilo:self.ilo+self.ng]

elif self.bc == "outflow":

# left boundary
self.u[0:self.ilo] = self.u[self.ilo]

# right boundary
self.u[self.ihi+1:] = self.u[self.ihi]

else:
sys.exit("invalid BC")


def norm(self, e):
""" return the norm of quantity e which lives on the grid """
if len(e) != 2*self.ng + self.nx:
return None

return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2))


class Simulation(object):

def __init__(self, grid):
self.grid = grid
self.t = 0.0


def init_cond(self, type="tophat"):

if type == "tophat":
self.grid.u[np.logical_and(self.grid.x >= 0.333,
self.grid.x <= 0.666)] = 1.0

elif type == "sine":
self.grid.u[:] = 1.0

index = np.logical_and(self.grid.x >= 0.333,
self.grid.x <= 0.666)
self.grid.u[index] += \
0.5*np.sin(2.0*np.pi*(self.grid.x[index]-0.333)/0.333)

elif type == "rarefaction":
self.grid.u[:] = 1.0
self.grid.u[self.grid.x > 0.5] = 2.0


def timestep(self, C):
return C*self.grid.dx/max(abs(self.grid.u[self.grid.ilo:
self.grid.ihi+1]))

def states(self, dt):
""" compute the left and right interface states """

g = self.grid
# compute the piecewise linear slopes -- 2nd order MC limiter
# we pick a range of cells that includes 1 ghost cell on either
# side
ib = g.ilo-1
ie = g.ihi+1

u = g.u

# this is the MC limiter from van Leer (1977), as given in
# LeVeque (2002). Note that this is slightly different than
# the expression from Colella (1990)

dc = g.scratch_array()
dl = g.scratch_array()
dr = g.scratch_array()

dc[ib:ie+1] = 0.5*(u[ib+1:ie+2] - u[ib-1:ie ])
dl[ib:ie+1] = u[ib+1:ie+2] - u[ib :ie+1]
dr[ib:ie+1] = u[ib :ie+1] - u[ib-1:ie ]

# these where's do a minmod()
d1 = 2.0*np.where(np.fabs(dl) < np.fabs(dr), dl, dr)
d2 = np.where(np.fabs(dc) < np.fabs(d1), dc, d1)
ldeltau = np.where(dl*dr > 0.0, d2, 0.0)

# now the interface states. Note that there are 1 more interfaces
# than zones
ul = g.scratch_array()
ur = g.scratch_array()

# are these indices right?
#
# --+-----------------+------------------+
# ^ i ^ ^ i+1
# ur(i) ul(i+1) ur(i+1)
#
ur[ib:ie+2] = u[ib:ie+2] - \
0.5*(1.0 + u[ib:ie+2]*dt/self.grid.dx)*ldeltau[ib:ie+2]

ul[ib+1:ie+2] = u[ib:ie+1] + \
0.5*(1.0 - u[ib:ie+1]*dt/self.grid.dx)*ldeltau[ib:ie+1]

return ul, ur


def riemann(self, ul, ur):
"""
Riemann problem for Burgers' equation.
"""

S = 0.5*(ul + ur)
ushock = np.where(S > 0.0, ul, ur)
ushock = np.where(S == 0.0, 0.0, ushock)

# rarefaction solution
urare = np.where(ur <= 0.0, ur, 0.0)
urare = np.where(ul >= 0.0, ul, urare)

us = np.where(ul > ur, ushock, urare)

return 0.5*us*us

def update(self, dt, flux):
""" conservative update """

g = self.grid

unew = g.scratch_array()

unew[g.ilo:g.ihi+1] = g.u[g.ilo:g.ihi+1] + \
dt/g.dx * (flux[g.ilo:g.ihi+1] - flux[g.ilo+1:g.ihi+2])

return unew

def evolve(self, C, tmax):

self.t = 0.0

g = self.grid

# main evolution loop
while (self.t < tmax):

# fill the boundary conditions
g.fill_BCs()

# get the timestep
dt = self.timestep(C)

if (self.t + dt > tmax):
dt = tmax - self.t

# get the interface states
ul, ur = self.states(dt)

# solve the Riemann problem at all interfaces
flux = self.riemann(ul, ur)

# do the conservative update
unew = self.update(dt, flux)

self.grid.u[:] = unew[:]

self.t += dt

def display_top(snapshot, key_type='lineno', limit=10):
snapshot = snapshot.filter_traces((
tracemalloc.Filter(False, "<frozen importlib._bootstrap>"),
tracemalloc.Filter(False, "<unknown>"),
))
top_stats = snapshot.statistics(key_type)

other = top_stats[limit:]
if other:
size = sum(stat.size for stat in other)
print("%s other: %.1f KiB" % (len(other), size / 1024))
total = sum(stat.size for stat in top_stats)
print("Total allocated size: %.1f KiB" % (total / 1024))

def main():
#-----------------------------------------------------------------------------
# sine

tracemalloc.start()

xmin = 0.0
xmax = 1.0
nx = 256
ng = 2
g = Grid1d(nx, ng, bc="periodic")

# maximum evolution time based on period for unit velocity
tmax = (xmax - xmin)/1.0

C = 0.8

plt.clf()

s = Simulation(g)

start = time.perf_counter()

for i in range(0, 10):
tend = (i+1)*0.02*tmax
s.init_cond("sine")

uinit = s.grid.u.copy()

s.evolve(C, tend)

c = 1.0 - (0.1 + i*0.1)
g = s.grid
#plt.plot(g.x[g.ilo:g.ihi+1], g.u[g.ilo:g.ihi+1], color=str(c))

timeTaken = time.perf_counter() - start

g = s.grid

plt.plot(g.x[g.ilo:g.ihi+1], uinit[g.ilo:g.ihi+1], ls=":", color="0.9", zorder=-1)

xs = g.x[g.ilo:g.ihi+1]
#us = uinit[g.ilo:g.ihi+1]
us = g.u[g.ilo:g.ihi+1]

plt.xlabel("$x$")
plt.ylabel("$u$")
plt.savefig("fv-burger-sine.pdf")

f = open('python-xs.txt','w')
for i in range(len(xs)):
f.write(str(xs[i]) + '\n')
f.close()

f = open('python-us.txt','w')
for i in range(len(us)):
f.write(str(us[i]) + '\n')
f.close()

snapshot = tracemalloc.take_snapshot()

display_top(snapshot)

print(asizeof.asizeof(g))
print(asizeof.asizeof(s))

tracemalloc.stop()

return timeTaken


#-----------------------------------------------------------------------------
# rarefaction
"""
xmin = 0.0
xmax = 1.0
nx = 256
ng = 2
g = Grid1d(nx, ng, bc="outflow")

# maximum evolution time based on period for unit velocity
tmax = (xmax - xmin)/1.0

C = 0.8

plt.clf()

s = Simulation(g)

for i in range(0, 10):
tend = (i+1)*0.02*tmax

s.init_cond("rarefaction")

uinit = s.grid.u.copy()

s.evolve(C, tend)

c = 1.0 - (0.1 + i*0.1)
plt.plot(g.x[g.ilo:g.ihi+1], g.u[g.ilo:g.ihi+1], color=str(c))


plt.plot(g.x[g.ilo:g.ihi+1], uinit[g.ilo:g.ihi+1], ls=":", color="0.9", zorder=-1)

plt.xlabel("$x$")
plt.ylabel("$u$")

plt.savefig("fv-burger-rarefaction.pdf")"""

if __name__ == "__main__":
main()

Binary file added Python Burgers example/fv-burger-sine.pdf
Binary file not shown.
10 changes: 10 additions & 0 deletions Python Burgers example/mprof0
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
Filename: C:\Users\alice\AppData\Local\Programs\Python\Python310\lib\site-packages\memory_profiler.py

Line # Mem usage Increment Occurrences Line Contents
=============================================================
1185 132.2 MiB 132.2 MiB 1 @wraps(wrapped=func)
1186 def wrapper(*args, **kwargs):
1187 132.2 MiB 0.0 MiB 1 prof = get_prof()
1188 132.2 MiB -0.0 MiB 1 val = prof(func)(*args, **kwargs)
1189 132.2 MiB 0.0 MiB 1 show_results_bound(prof)
1190 132.2 MiB 0.0 MiB 1 return val
4 changes: 4 additions & 0 deletions Python Burgers example/mprofile_20230427223539.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
CMDLINE C:\Users\alice\AppData\Local\Programs\Python\Python310\python.exe run-burgers.py
MEM 1.593750 1682631339.4155
MEM 12.250000 1682631339.5263
MEM 18.988281 1682631339.6373
Loading