From 9501363c74f5b1893a46fd020cddc982f66a4a29 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Mon, 11 Feb 2019 11:59:56 +0000 Subject: [PATCH 01/33] Start on DG methods --- advection/dg.py | 216 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 advection/dg.py diff --git a/advection/dg.py b/advection/dg.py new file mode 100644 index 0000000..3d66ae0 --- /dev/null +++ b/advection/dg.py @@ -0,0 +1,216 @@ +""" +Discontinuous Galerkin for the advection equation. +""" + +import numpy +from matplotlib import pyplot +import matplotlib as mpl +import quadpy + +mpl.rcParams['mathtext.fontset'] = 'cm' +mpl.rcParams['mathtext.rm'] = 'serif' + + + +class Grid1d(object): + + def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): + + self.ng = ng + self.nx = nx + self.np = np + + self.xmin = xmin + self.xmax = xmax + + # 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 + (numpy.arange(nx+2*ng)-ng+0.5)*self.dx + self.xl = xmin + (numpy.arange(nx+2*ng)-ng)*self.dx + self.xr = xmin + (numpy.arange(nx+2*ng)+1.0)*self.dx + + # storage for the solution + # These are the modes of the solution at each point, so the + # first index is the mode + self.a = numpy.zeros((self.np, nx+2*ng), dtype=numpy.float64) + + # Need the Gauss-Lobatto nodes and weights in the reference element + GL = quadpy.line_segment.GaussLobatto(np+2) + self.nodes = GL.points + self.weights = GL.weights + # To go from modal to nodal we need the Vandermonde matrix + self.V = numpy.zeros((np+2, np+2)) + for n in range(np+2): + c = numpy.zeros(n) + c[n] = 1 + self.V[:, n] = numpy.polynomial.legendre.legval(self.nodes, c) + # We can now get the nodal values from self.V @ self.a[:, i] + + # Need the weights multiplied by P_p' for the interior flux + self.modified_weights = numpy.zeros((np, np+2)) + for p in range(np): + pp_c = numpy.zeros(p+1) + pp_c[p] = 1 + pp_prime_c = numpy.polynomial.legendre.legder(pp_c) + pp_prime_nodes = numpy.polynomial.legendre.legval(self.nodes, pp_prime_c) + self.modified_weights[p, :] = self.weights * pp_prime_nodes + + # Nodes in the computational coordinates + self.all_nodes = numpy.zeros((np+2)*(nx+2*ng), dtype=numpy.float64) + for i in range(nx+2*ng): + self.all_nodes[(np+2)*i: (np+2)*(i+1)] = self.x + self.nodes * self.dx / 2 + + def scratch_array(self): + """ return a scratch array dimensioned for our grid """ + return numpy.zeros((self.np, self.nx+2*self.ng), dtype=numpy.float64) + + + def fill_BCs(self): + """ fill all single ghostcell with periodic boundary conditions """ + + for n in range(self.ng): + # left boundary + self.a[:, self.ilo-1-n] = self.a[:, self.ihi-n] + + # right boundary + self.a[:, self.ihi+1+n] = self.a[:, self.ilo+n] + + 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 numpy.sqrt(self.dx*numpy.sum(e[self.ilo:self.ihi+1]**2)) + return numpy.max(abs(e[0, self.ilo:self.ihi+1])) + + +class Simulation(object): + + def __init__(self, grid, u, C=0.8): + self.grid = grid + self.t = 0.0 # simulation time + self.u = u # the constant advective velocity + self.C = C # CFL number + + + def init_cond(self, type="sine"): + """ initialize the data """ + if type == "tophat": + init_a = lambda x : numpy.where(numpy.logical_and(x >=0.333, x <=0.666), + numpy.ones_like(x), numpy.zeros_like(x)) + + elif type == "sine": + init_a = lambda x : numpy.sin(2.0*numpy.pi*x/(self.grid.xmax-self.grid.xmin)) + + nodal_a = init_a(self.grid.all_nodes) + for i in range(self.grid.np+2*self.grid.ng): + for p in range(self.grid.np): + self.grid.a[p, i] = numpy.sum(nodal_a[(self.grid.np+2)*i:(self.grid.np+2)*(i+1)] * self.V[:,p]) + + + def timestep(self): + """ return the advective timestep """ + return self.C*self.grid.dx/self.u + + + def period(self): + """ return the period for advection with velocity u """ + return (self.grid.xmax - self.grid.xmin)/self.u + + + def states(self): + """ compute the left and right interface states """ + + # Evaluate the nodal values at the domain edges + g = self.grid + + al = g.scratch_array() + ar = g.scratch_array() + + # i is looping over interfaces, so al is the right edge of the left + # element, etc. + for i in range(g.ilo, g.ihi+2): + for p in range(g.np): + al[p, i] = numpy.dot(g.a[:, i-1], g.V[-1, :]) + ar[p, i] = numpy.dot(g.a[:, i ], g.V[ 0, :]) + + return al, ar + + + def riemann(self, al, ar): + """ + Riemann problem for advection -- this is simply upwinding, + but we return the flux + """ + + if self.u > 0.0: + return self.u*al + else: + return self.u*ar + + + def rk_substep(self): + + g = self.grid + g.fill_BCs() + rhs = g.scratch_array() + + # Integrate flux over element + interior_f = g.scratch_array() + for p in range(g.np): + for i in range(g.ilo, g.ihi+1): + nodal_a = g.V @ g.a[:, i] + nodal_f = self.u * nodal_a + interior_f[p, i] = numpy.dot(nodal_f, g.modified_weights[p, :]) + + # Use Riemann solver to get fluxes between elements + boundary_f = self.riemann(*self.states()) + rhs = interior_f + for p in range(g.np): + for i in range(g.ilo, g.ihi+1): + rhs[p, i] += (-1)**p * boundary_f[p, i] - boundary_f[p, i+1] + + # Multiply by mass matrix, which is diagonal. + for p in range(g.np): + rhs[p, :] *= (2*p + 1) / 2 + + return rhs + + + def evolve(self, num_periods=1): + """ evolve the linear advection equation using RK3 (SSP) """ + self.t = 0.0 + g = self.grid + + tmax = num_periods*self.period() + + # main evolution loop + while self.t < tmax: + + # fill the boundary conditions + g.fill_BCs() + + # get the timestep + dt = self.timestep() + + if self.t + dt > tmax: + dt = tmax - self.t + + # RK3 - SSP + # Store the data at the start of the step + a_start = g.a.copy() + k1 = dt * self.rk_substep() + g.a = a_start + k1 + # a1 = g.a.copy() + k2 = dt * self.rk_substep() + g.a = a_start + k2 / 4 + a2 = g.a.copy() + k3 = dt * self.rk_substep() + g.a = (a_start + 2 * a2 + 2 * k3) / 3 + + self.t += dt From 924457b7ad5ff167074e339f809432c29a0e552f Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Mon, 11 Feb 2019 16:54:37 +0000 Subject: [PATCH 02/33] Fix plotting issue, nodal to modal --- advection/dg.py | 76 ++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 62 insertions(+), 14 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 3d66ae0..38565b1 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -16,6 +16,8 @@ class Grid1d(object): def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): + assert np > 1 + self.ng = ng self.nx = nx self.np = np @@ -40,19 +42,24 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): self.a = numpy.zeros((self.np, nx+2*ng), dtype=numpy.float64) # Need the Gauss-Lobatto nodes and weights in the reference element - GL = quadpy.line_segment.GaussLobatto(np+2) + GL = quadpy.line_segment.GaussLobatto(np) self.nodes = GL.points self.weights = GL.weights # To go from modal to nodal we need the Vandermonde matrix - self.V = numpy.zeros((np+2, np+2)) - for n in range(np+2): - c = numpy.zeros(n) - c[n] = 1 - self.V[:, n] = numpy.polynomial.legendre.legval(self.nodes, c) - # We can now get the nodal values from self.V @ self.a[:, i] + self.V = numpy.zeros((np, np)) + c = numpy.eye(np) + self.V = numpy.polynomial.legendre.legval(self.nodes, c) + self.V_inv = numpy.linalg.inv(self.V) +# print("first", self.V) +# for n in range(np): +# c = numpy.zeros(n+1) +# c[n] = 1 +# self.V[:, n] = numpy.polynomial.legendre.legval(self.nodes, c) +# # We can now get the nodal values from self.V @ self.a[:, i] +# print("second", self.V) # Need the weights multiplied by P_p' for the interior flux - self.modified_weights = numpy.zeros((np, np+2)) + self.modified_weights = numpy.zeros((np, np)) for p in range(np): pp_c = numpy.zeros(p+1) pp_c[p] = 1 @@ -61,9 +68,25 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): self.modified_weights[p, :] = self.weights * pp_prime_nodes # Nodes in the computational coordinates - self.all_nodes = numpy.zeros((np+2)*(nx+2*ng), dtype=numpy.float64) + self.all_nodes = numpy.zeros((np)*(nx+2*ng), dtype=numpy.float64) + self.all_nodes_per_node = numpy.zeros_like(self.a) for i in range(nx+2*ng): - self.all_nodes[(np+2)*i: (np+2)*(i+1)] = self.x + self.nodes * self.dx / 2 + self.all_nodes[(np)*i: (np)*(i+1)] = self.x[i] + self.nodes * self.dx / 2 + self.all_nodes_per_node[:, i] = self.x[i] + self.nodes * self.dx / 2 + + def modal_to_nodal(self): + nodal = numpy.zeros_like(self.a) + for i in range(self.nx+2*self.ng): + nodal[:, i] = self.V @ self.a[:, i] + return nodal + + def nodal_to_modal(self, nodal): + for i in range(self.nx+2*self.ng): + self.a[:, i] = self.V_inv @ nodal[:, i] + + def plotting_data(self): + return (self.all_nodes, + self.modal_to_nodal().ravel(order='F')) def scratch_array(self): """ return a scratch array dimensioned for our grid """ @@ -107,10 +130,12 @@ def init_cond(self, type="sine"): elif type == "sine": init_a = lambda x : numpy.sin(2.0*numpy.pi*x/(self.grid.xmax-self.grid.xmin)) - nodal_a = init_a(self.grid.all_nodes) - for i in range(self.grid.np+2*self.grid.ng): - for p in range(self.grid.np): - self.grid.a[p, i] = numpy.sum(nodal_a[(self.grid.np+2)*i:(self.grid.np+2)*(i+1)] * self.V[:,p]) + nodal_a = init_a(self.grid.all_nodes_per_node) +# pyplot.plot(self.grid.all_nodes, nodal_a.ravel(order='F')) +# pyplot.show() + self.grid.nodal_to_modal(nodal_a) +# for i in range(self.grid.nx+2*self.grid.ng): +# self.grid.a[:, i] = numpy.linalg.inv(self.grid.V) @ nodal_a[(self.grid.np)*i:(self.grid.np)*(i+1)] def timestep(self): @@ -214,3 +239,26 @@ def evolve(self, num_periods=1): g.a = (a_start + 2 * a2 + 2 * k3) / 3 self.t += dt + + +if __name__ == "__main__": + + + #------------------------------------------------------------------------- + # compare limiting and no-limiting + + xmin = 0.0 + xmax = 1.0 + nx = 500 + ng = 1 + + g = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=5) + + u = 1.0 + + s = Simulation(g, u, C=0.5) + s.init_cond("sine") + ainit = s.grid.a.copy() + plot_x, plot_a = g.plotting_data() + + pyplot.plot(plot_x, plot_a, 'kx') \ No newline at end of file From a8c7b0c0fac9dac9d35abef8e91fc072b2568b7f Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 12 Feb 2019 11:09:02 +0000 Subject: [PATCH 03/33] Got to the evolution, but check mass matrix normalization --- advection/dg.py | 149 +++++++++++++++++++++++++++--------------------- 1 file changed, 85 insertions(+), 64 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 38565b1..a4b2939 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -11,13 +11,12 @@ mpl.rcParams['mathtext.rm'] = 'serif' - class Grid1d(object): def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): assert np > 1 - + self.ng = ng self.nx = nx self.np = np @@ -40,7 +39,7 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): # These are the modes of the solution at each point, so the # first index is the mode self.a = numpy.zeros((self.np, nx+2*ng), dtype=numpy.float64) - + # Need the Gauss-Lobatto nodes and weights in the reference element GL = quadpy.line_segment.GaussLobatto(np) self.nodes = GL.points @@ -50,40 +49,36 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): c = numpy.eye(np) self.V = numpy.polynomial.legendre.legval(self.nodes, c) self.V_inv = numpy.linalg.inv(self.V) -# print("first", self.V) -# for n in range(np): -# c = numpy.zeros(n+1) -# c[n] = 1 -# self.V[:, n] = numpy.polynomial.legendre.legval(self.nodes, c) -# # We can now get the nodal values from self.V @ self.a[:, i] -# print("second", self.V) - + # Need the weights multiplied by P_p' for the interior flux self.modified_weights = numpy.zeros((np, np)) for p in range(np): pp_c = numpy.zeros(p+1) pp_c[p] = 1 pp_prime_c = numpy.polynomial.legendre.legder(pp_c) - pp_prime_nodes = numpy.polynomial.legendre.legval(self.nodes, pp_prime_c) + pp_prime_nodes = numpy.polynomial.legendre.legval(self.nodes, + pp_prime_c) self.modified_weights[p, :] = self.weights * pp_prime_nodes - + # Nodes in the computational coordinates self.all_nodes = numpy.zeros((np)*(nx+2*ng), dtype=numpy.float64) self.all_nodes_per_node = numpy.zeros_like(self.a) for i in range(nx+2*ng): - self.all_nodes[(np)*i: (np)*(i+1)] = self.x[i] + self.nodes * self.dx / 2 - self.all_nodes_per_node[:, i] = self.x[i] + self.nodes * self.dx / 2 - + self.all_nodes[(np)*i:(np)*(i+1)] = (self.x[i] + + self.nodes * self.dx / 2) + self.all_nodes_per_node[:, i] = (self.x[i] + + self.nodes * self.dx / 2) + def modal_to_nodal(self): nodal = numpy.zeros_like(self.a) for i in range(self.nx+2*self.ng): nodal[:, i] = self.V @ self.a[:, i] return nodal - + def nodal_to_modal(self, nodal): for i in range(self.nx+2*self.ng): self.a[:, i] = self.V_inv @ nodal[:, i] - + def plotting_data(self): return (self.all_nodes, self.modal_to_nodal().ravel(order='F')) @@ -92,14 +87,12 @@ def scratch_array(self): """ return a scratch array dimensioned for our grid """ return numpy.zeros((self.np, self.nx+2*self.ng), dtype=numpy.float64) - def fill_BCs(self): """ fill all single ghostcell with periodic boundary conditions """ for n in range(self.ng): # left boundary self.a[:, self.ilo-1-n] = self.a[:, self.ihi-n] - # right boundary self.a[:, self.ihi+1+n] = self.a[:, self.ilo+n] @@ -108,7 +101,7 @@ def norm(self, e): if len(e) != 2*self.ng + self.nx: return None - #return numpy.sqrt(self.dx*numpy.sum(e[self.ilo:self.ihi+1]**2)) + # return numpy.sqrt(self.dx*numpy.sum(e[self.ilo:self.ihi+1]**2)) return numpy.max(abs(e[0, self.ilo:self.ihi+1])) @@ -116,57 +109,54 @@ class Simulation(object): def __init__(self, grid, u, C=0.8): self.grid = grid - self.t = 0.0 # simulation time - self.u = u # the constant advective velocity - self.C = C # CFL number - + self.t = 0.0 # simulation time + self.u = u # the constant advective velocity + self.C = C # CFL number def init_cond(self, type="sine"): """ initialize the data """ if type == "tophat": - init_a = lambda x : numpy.where(numpy.logical_and(x >=0.333, x <=0.666), - numpy.ones_like(x), numpy.zeros_like(x)) + def init_a(x): + return numpy.where(numpy.logical_and(x >= 0.333, + x <= 0.666), + numpy.ones_like(x), + numpy.zeros_like(x)) elif type == "sine": - init_a = lambda x : numpy.sin(2.0*numpy.pi*x/(self.grid.xmax-self.grid.xmin)) + def init_a(x): + return numpy.sin(2.0 * numpy.pi * x / + (self.grid.xmax - self.grid.xmin)) nodal_a = init_a(self.grid.all_nodes_per_node) -# pyplot.plot(self.grid.all_nodes, nodal_a.ravel(order='F')) -# pyplot.show() self.grid.nodal_to_modal(nodal_a) -# for i in range(self.grid.nx+2*self.grid.ng): -# self.grid.a[:, i] = numpy.linalg.inv(self.grid.V) @ nodal_a[(self.grid.np)*i:(self.grid.np)*(i+1)] - def timestep(self): """ return the advective timestep """ return self.C*self.grid.dx/self.u - def period(self): """ return the period for advection with velocity u """ return (self.grid.xmax - self.grid.xmin)/self.u - def states(self): """ compute the left and right interface states """ # Evaluate the nodal values at the domain edges g = self.grid - al = g.scratch_array() - ar = g.scratch_array() + al = numpy.zeros(g.nx+2*g.ng) + ar = numpy.zeros(g.nx+2*g.ng) + + nodal = g.modal_to_nodal() # i is looping over interfaces, so al is the right edge of the left # element, etc. for i in range(g.ilo, g.ihi+2): - for p in range(g.np): - al[p, i] = numpy.dot(g.a[:, i-1], g.V[-1, :]) - ar[p, i] = numpy.dot(g.a[:, i ], g.V[ 0, :]) + al[i] = nodal[-1, i-1] + ar[i] = nodal[ 0, i ] return al, ar - def riemann(self, al, ar): """ Riemann problem for advection -- this is simply upwinding, @@ -178,13 +168,14 @@ def riemann(self, al, ar): else: return self.u*ar - def rk_substep(self): - + """ + Take a single RK substep + """ g = self.grid g.fill_BCs() rhs = g.scratch_array() - + # Integrate flux over element interior_f = g.scratch_array() for p in range(g.np): @@ -192,20 +183,19 @@ def rk_substep(self): nodal_a = g.V @ g.a[:, i] nodal_f = self.u * nodal_a interior_f[p, i] = numpy.dot(nodal_f, g.modified_weights[p, :]) - # Use Riemann solver to get fluxes between elements boundary_f = self.riemann(*self.states()) rhs = interior_f for p in range(g.np): for i in range(g.ilo, g.ihi+1): - rhs[p, i] += (-1)**p * boundary_f[p, i] - boundary_f[p, i+1] - - # Multiply by mass matrix, which is diagonal. - for p in range(g.np): - rhs[p, :] *= (2*p + 1) / 2 - - return rhs + rhs[p, i] += (-1)**p * boundary_f[i] - boundary_f[i+1] + + # Multiply by mass matrix (inverse), which is diagonal. + # Is it orthonormal? +# for p in range(g.np): +# rhs[p, :] *= (2*p + 1) / 2 + return rhs def evolve(self, num_periods=1): """ evolve the linear advection equation using RK3 (SSP) """ @@ -216,7 +206,6 @@ def evolve(self, num_periods=1): # main evolution loop while self.t < tmax: - # fill the boundary conditions g.fill_BCs() @@ -243,22 +232,54 @@ def evolve(self, num_periods=1): if __name__ == "__main__": - - #------------------------------------------------------------------------- - # compare limiting and no-limiting + # ------------------------------------------------------------------------- + # DG of sine wave xmin = 0.0 xmax = 1.0 - nx = 500 + nx = 4 ng = 1 - g = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=5) + g3 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=3) + g7 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=7) u = 1.0 - s = Simulation(g, u, C=0.5) - s.init_cond("sine") - ainit = s.grid.a.copy() - plot_x, plot_a = g.plotting_data() - - pyplot.plot(plot_x, plot_a, 'kx') \ No newline at end of file + # The CFL limit for DG is reduced by a factor 1/(2 p + 1) + s3 = Simulation(g3, u, C=0.8/(2*3+1)) + s3.init_cond("sine") + s7 = Simulation(g7, u, C=0.8/(2*7+1)) + s7.init_cond("sine") + # Plot the initial data to show how, difference in nodal locations as + # number of modes varies + plot_x3, plot_a3 = g3.plotting_data() + a3init = plot_a3.copy() + plot_x7, plot_a7 = g7.plotting_data() + a7init = plot_a7.copy() + pyplot.plot(plot_x3, plot_a3, 'bo', label=r"$p=3$") + pyplot.plot(plot_x7, plot_a7, 'r^', label=r"$p=7$") + pyplot.xlim(0, 1) + pyplot.xlabel(r'$x$') + pyplot.ylabel(r'$a$') + pyplot.legend() + pyplot.show() + + s3.evolve(num_periods=1) + plot_x3, plot_a3 = g3.plotting_data() + pyplot.plot(plot_x3, plot_a3, 'bo', label=r"$p=3$") + pyplot.plot(plot_x3, a3init, 'r^', label=r"$p=3$, initial") + pyplot.xlim(0, 1) + pyplot.xlabel(r'$x$') + pyplot.ylabel(r'$a$') + pyplot.legend() + pyplot.show() + + s7.evolve(num_periods=1) + plot_x7, plot_a7 = g7.plotting_data() + pyplot.plot(plot_x7, plot_a7, 'bo', label=r"$p=7$") + pyplot.plot(plot_x7, a7init, 'r^', label=r"$p=7$, initial") + pyplot.xlim(0, 1) + pyplot.xlabel(r'$x$') + pyplot.ylabel(r'$a$') + pyplot.legend() + pyplot.show() From 7130a5097ad3719c792abd1434dbee73fe911a3b Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 12 Feb 2019 11:34:18 +0000 Subject: [PATCH 04/33] A different mass matrix, which isn't diagonal? --- advection/dg.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/advection/dg.py b/advection/dg.py index a4b2939..e2cedb9 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -49,6 +49,8 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): c = numpy.eye(np) self.V = numpy.polynomial.legendre.legval(self.nodes, c) self.V_inv = numpy.linalg.inv(self.V) + self.M_inv = 2 / self.dx * self.V @ self.V.T + print("Inverse mass matrix is", self.M_inv) # Need the weights multiplied by P_p' for the interior flux self.modified_weights = numpy.zeros((np, np)) @@ -194,6 +196,9 @@ def rk_substep(self): # Is it orthonormal? # for p in range(g.np): # rhs[p, :] *= (2*p + 1) / 2 +# for p in range(g.np): +# rhs[p, :] *= (2 * p + 1) / g.dx**p + rhs = g.M_inv @ rhs return rhs From 436ff4fa284030def029781025ccee705e2fce09 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Thu, 14 Feb 2019 13:16:48 +0000 Subject: [PATCH 05/33] Now works for linear smooth wave --- advection/dg.py | 103 +++++++++++++++++++++++++----------------------- 1 file changed, 53 insertions(+), 50 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index e2cedb9..99fa8e7 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -15,7 +15,7 @@ class Grid1d(object): def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): - assert np > 1 + assert np > 0 self.ng = ng self.nx = nx @@ -38,36 +38,36 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): # storage for the solution # These are the modes of the solution at each point, so the # first index is the mode - self.a = numpy.zeros((self.np, nx+2*ng), dtype=numpy.float64) + # NO! These are actually the *nodal* values + self.a = numpy.zeros((self.np+1, nx+2*ng), dtype=numpy.float64) # Need the Gauss-Lobatto nodes and weights in the reference element - GL = quadpy.line_segment.GaussLobatto(np) + GL = quadpy.line_segment.GaussLobatto(np+1) self.nodes = GL.points self.weights = GL.weights # To go from modal to nodal we need the Vandermonde matrix - self.V = numpy.zeros((np, np)) - c = numpy.eye(np) - self.V = numpy.polynomial.legendre.legval(self.nodes, c) + self.V = numpy.polynomial.legendre.legvander(self.nodes, np) + c = numpy.eye(np+1) + # Orthonormalize + for p in range(np+1): + self.V[:, p] /= numpy.sqrt(2/(2*p+1)) + c[p, p] /= numpy.sqrt(2/(2*p+1)) self.V_inv = numpy.linalg.inv(self.V) - self.M_inv = 2 / self.dx * self.V @ self.V.T - print("Inverse mass matrix is", self.M_inv) - - # Need the weights multiplied by P_p' for the interior flux - self.modified_weights = numpy.zeros((np, np)) - for p in range(np): - pp_c = numpy.zeros(p+1) - pp_c[p] = 1 - pp_prime_c = numpy.polynomial.legendre.legder(pp_c) - pp_prime_nodes = numpy.polynomial.legendre.legval(self.nodes, - pp_prime_c) - self.modified_weights[p, :] = self.weights * pp_prime_nodes + self.M = numpy.linalg.inv(self.V @ self.V.T) + self.M_inv = self.V @ self.V.T + # Derivatives of Legendre polynomials lead to derivatives of V + dV = numpy.polynomial.legendre.legval(self.nodes, + numpy.polynomial.legendre.legder(c)).T + self.D = dV @ self.V_inv + # Stiffness matrix for the interior flux + self.S = self.M @ self.D # Nodes in the computational coordinates - self.all_nodes = numpy.zeros((np)*(nx+2*ng), dtype=numpy.float64) + self.all_nodes = numpy.zeros((np+1)*(nx+2*ng), dtype=numpy.float64) self.all_nodes_per_node = numpy.zeros_like(self.a) for i in range(nx+2*ng): - self.all_nodes[(np)*i:(np)*(i+1)] = (self.x[i] + - self.nodes * self.dx / 2) + self.all_nodes[(np+1)*i:(np+1)*(i+1)] = (self.x[i] + + self.nodes * self.dx / 2) self.all_nodes_per_node[:, i] = (self.x[i] + self.nodes * self.dx / 2) @@ -83,11 +83,11 @@ def nodal_to_modal(self, nodal): def plotting_data(self): return (self.all_nodes, - self.modal_to_nodal().ravel(order='F')) + self.a.ravel(order='F')) def scratch_array(self): """ return a scratch array dimensioned for our grid """ - return numpy.zeros((self.np, self.nx+2*self.ng), dtype=numpy.float64) + return numpy.zeros((self.np+1, self.nx+2*self.ng), dtype=numpy.float64) def fill_BCs(self): """ fill all single ghostcell with periodic boundary conditions """ @@ -129,8 +129,7 @@ def init_a(x): return numpy.sin(2.0 * numpy.pi * x / (self.grid.xmax - self.grid.xmin)) - nodal_a = init_a(self.grid.all_nodes_per_node) - self.grid.nodal_to_modal(nodal_a) + self.grid.a = init_a(self.grid.all_nodes_per_node) def timestep(self): """ return the advective timestep """ @@ -149,13 +148,11 @@ def states(self): al = numpy.zeros(g.nx+2*g.ng) ar = numpy.zeros(g.nx+2*g.ng) - nodal = g.modal_to_nodal() - # i is looping over interfaces, so al is the right edge of the left # element, etc. for i in range(g.ilo, g.ihi+2): - al[i] = nodal[-1, i-1] - ar[i] = nodal[ 0, i ] + al[i] = g.a[-1, i-1] + ar[i] = g.a[ 0, i ] return al, ar @@ -179,28 +176,18 @@ def rk_substep(self): rhs = g.scratch_array() # Integrate flux over element - interior_f = g.scratch_array() - for p in range(g.np): - for i in range(g.ilo, g.ihi+1): - nodal_a = g.V @ g.a[:, i] - nodal_f = self.u * nodal_a - interior_f[p, i] = numpy.dot(nodal_f, g.modified_weights[p, :]) + f = self.u * g.a + interior_f = g.S.T @ f # Use Riemann solver to get fluxes between elements boundary_f = self.riemann(*self.states()) rhs = interior_f - for p in range(g.np): - for i in range(g.ilo, g.ihi+1): - rhs[p, i] += (-1)**p * boundary_f[i] - boundary_f[i+1] + rhs[ 0, 1:-1] += boundary_f[1:-1] + rhs[-1, 1:-1] -= boundary_f[2:] - # Multiply by mass matrix (inverse), which is diagonal. - # Is it orthonormal? -# for p in range(g.np): -# rhs[p, :] *= (2*p + 1) / 2 -# for p in range(g.np): -# rhs[p, :] *= (2 * p + 1) / g.dx**p - rhs = g.M_inv @ rhs + # Multiply by mass matrix (inverse). + rhs_i = 2 / g.dx * g.M_inv @ rhs - return rhs + return rhs_i def evolve(self, num_periods=1): """ evolve the linear advection equation using RK3 (SSP) """ @@ -225,9 +212,9 @@ def evolve(self, num_periods=1): a_start = g.a.copy() k1 = dt * self.rk_substep() g.a = a_start + k1 - # a1 = g.a.copy() + a1 = g.a.copy() k2 = dt * self.rk_substep() - g.a = a_start + k2 / 4 + g.a = (3 * a_start + a1 + k2) / 4 a2 = g.a.copy() k3 = dt * self.rk_substep() g.a = (a_start + 2 * a2 + 2 * k3) / 3 @@ -242,25 +229,31 @@ def evolve(self, num_periods=1): xmin = 0.0 xmax = 1.0 - nx = 4 + nx = 16 ng = 1 + g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=1) g3 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=3) g7 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=7) u = 1.0 # The CFL limit for DG is reduced by a factor 1/(2 p + 1) + s1 = Simulation(g1, u, C=0.8/(2*1+1)) + s1.init_cond("sine") s3 = Simulation(g3, u, C=0.8/(2*3+1)) s3.init_cond("sine") - s7 = Simulation(g7, u, C=0.8/(2*7+1)) + s7 = Simulation(g7, u, C=0.1/(2*7+1)) # Not sure what the critical CFL is s7.init_cond("sine") # Plot the initial data to show how, difference in nodal locations as # number of modes varies + plot_x1, plot_a1 = g1.plotting_data() + a1init = plot_a1.copy() plot_x3, plot_a3 = g3.plotting_data() a3init = plot_a3.copy() plot_x7, plot_a7 = g7.plotting_data() a7init = plot_a7.copy() + pyplot.plot(plot_x1, plot_a1, 'k>', label=r"$p=1$") pyplot.plot(plot_x3, plot_a3, 'bo', label=r"$p=3$") pyplot.plot(plot_x7, plot_a7, 'r^', label=r"$p=7$") pyplot.xlim(0, 1) @@ -269,6 +262,16 @@ def evolve(self, num_periods=1): pyplot.legend() pyplot.show() + s1.evolve(num_periods=1) + plot_x1, plot_a1 = g1.plotting_data() + pyplot.plot(plot_x1, plot_a1, 'k>', label=r"$p=1$") + pyplot.plot(plot_x1, a1init, 'r^', label=r"$p=1$, initial") + pyplot.xlim(0, 1) + pyplot.xlabel(r'$x$') + pyplot.ylabel(r'$a$') + pyplot.legend() + pyplot.show() + s3.evolve(num_periods=1) plot_x3, plot_a3 = g3.plotting_data() pyplot.plot(plot_x3, plot_a3, 'bo', label=r"$p=3$") From a0339da76f612fd15c980cc7b7134e5c5b83c64f Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Thu, 14 Feb 2019 17:21:20 +0000 Subject: [PATCH 06/33] Start looking at limiting --- advection/dg.py | 50 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 37 insertions(+), 13 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 99fa8e7..35bb291 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -71,15 +71,17 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): self.all_nodes_per_node[:, i] = (self.x[i] + self.nodes * self.dx / 2) - def modal_to_nodal(self): - nodal = numpy.zeros_like(self.a) +# def modal_to_nodal(self): +# nodal = numpy.zeros_like(self.a) +# for i in range(self.nx+2*self.ng): +# nodal[:, i] = self.V @ self.a[:, i] +# return nodal + + def nodal_to_modal(self): + modal = numpy.zeros_like(self.a) for i in range(self.nx+2*self.ng): - nodal[:, i] = self.V @ self.a[:, i] - return nodal - - def nodal_to_modal(self, nodal): - for i in range(self.nx+2*self.ng): - self.a[:, i] = self.V_inv @ nodal[:, i] + modal[:, i] = self.V_inv @ self.a[:, i] + return modal def plotting_data(self): return (self.all_nodes, @@ -154,6 +156,28 @@ def states(self): al[i] = g.a[-1, i-1] ar[i] = g.a[ 0, i ] + # TODO: remove this return, fix the limiting. + return al, ar + + # Limiting, using minmod (Hesthaven p 443-4) + theta = 2 + a_modal = g.nodal_to_modal() + a_average = a_modal[0, :] + a_m = numpy.zeros_like(a_average) + a_m[1:] = a_average[:-1] + a_el = g.a[0 , :] + a_er = g.a[-1, :] + + check_left = a_average - minmod([a_average - a_el, + a_average - a_m, + a_er - a_average]) + check_right = a_average + minmod([a_average - a_el, + a_average - m, + a_er - a_average]) + ids = numpy.logical_or(numpy.isclose(a_el, check_left), + numpy.isclose(a_er, check_right)) + + return al, ar def riemann(self, al, ar): @@ -229,7 +253,7 @@ def evolve(self, num_periods=1): xmin = 0.0 xmax = 1.0 - nx = 16 + nx = 64 ng = 1 g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=1) @@ -240,11 +264,11 @@ def evolve(self, num_periods=1): # The CFL limit for DG is reduced by a factor 1/(2 p + 1) s1 = Simulation(g1, u, C=0.8/(2*1+1)) - s1.init_cond("sine") s3 = Simulation(g3, u, C=0.8/(2*3+1)) - s3.init_cond("sine") - s7 = Simulation(g7, u, C=0.1/(2*7+1)) # Not sure what the critical CFL is - s7.init_cond("sine") + s7 = Simulation(g7, u, C=0.5/(2*7+1)) # Not sure what the critical CFL is + for s in [s1, s3, s7]: +# s.init_cond("sine") + s.init_cond("tophat") # Plot the initial data to show how, difference in nodal locations as # number of modes varies plot_x1, plot_a1 = g1.plotting_data() From da086b20d37fda333046db85384bee4ed9f9ada3 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 26 Feb 2019 14:59:20 +0000 Subject: [PATCH 07/33] Convergence for low m --- advection/dg.py | 83 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 62 insertions(+), 21 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 35bb291..5cd0a52 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -13,13 +13,13 @@ class Grid1d(object): - def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): + def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): - assert np > 0 + assert m > 0 self.ng = ng self.nx = nx - self.np = np + self.m = m self.xmin = xmin self.xmax = xmax @@ -39,17 +39,17 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): # These are the modes of the solution at each point, so the # first index is the mode # NO! These are actually the *nodal* values - self.a = numpy.zeros((self.np+1, nx+2*ng), dtype=numpy.float64) + self.a = numpy.zeros((self.m+1, nx+2*ng), dtype=numpy.float64) # Need the Gauss-Lobatto nodes and weights in the reference element - GL = quadpy.line_segment.GaussLobatto(np+1) + GL = quadpy.line_segment.GaussLobatto(m+1) self.nodes = GL.points self.weights = GL.weights # To go from modal to nodal we need the Vandermonde matrix - self.V = numpy.polynomial.legendre.legvander(self.nodes, np) - c = numpy.eye(np+1) + self.V = numpy.polynomial.legendre.legvander(self.nodes, m) + c = numpy.eye(m+1) # Orthonormalize - for p in range(np+1): + for p in range(m+1): self.V[:, p] /= numpy.sqrt(2/(2*p+1)) c[p, p] /= numpy.sqrt(2/(2*p+1)) self.V_inv = numpy.linalg.inv(self.V) @@ -63,10 +63,10 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, np=3): self.S = self.M @ self.D # Nodes in the computational coordinates - self.all_nodes = numpy.zeros((np+1)*(nx+2*ng), dtype=numpy.float64) + self.all_nodes = numpy.zeros((m+1)*(nx+2*ng), dtype=numpy.float64) self.all_nodes_per_node = numpy.zeros_like(self.a) for i in range(nx+2*ng): - self.all_nodes[(np+1)*i:(np+1)*(i+1)] = (self.x[i] + + self.all_nodes[(m+1)*i:(m+1)*(i+1)] = (self.x[i] + self.nodes * self.dx / 2) self.all_nodes_per_node[:, i] = (self.x[i] + self.nodes * self.dx / 2) @@ -89,7 +89,7 @@ def plotting_data(self): def scratch_array(self): """ return a scratch array dimensioned for our grid """ - return numpy.zeros((self.np+1, self.nx+2*self.ng), dtype=numpy.float64) + return numpy.zeros((self.m+1, self.nx+2*self.ng), dtype=numpy.float64) def fill_BCs(self): """ fill all single ghostcell with periodic boundary conditions """ @@ -101,12 +101,24 @@ def fill_BCs(self): self.a[:, self.ihi+1+n] = self.a[:, self.ilo+n] def norm(self, e): - """ return the norm of quantity e which lives on the grid """ - if len(e) != 2*self.ng + self.nx: + """ + Return the norm of quantity e which lives on the grid. + + This is the 'broken norm': the quantity is integrated over each + individual element using Gauss-Lobatto quadrature (as we have those + nodes and weights), and the 2-norm of the result is then returned. + """ + if not numpy.allclose(e.shape, self.all_nodes_per_node.shape): return None - # return numpy.sqrt(self.dx*numpy.sum(e[self.ilo:self.ihi+1]**2)) - return numpy.max(abs(e[0, self.ilo:self.ihi+1])) + + # This is L_inf norm... +# return numpy.max(abs(e[:, self.ilo:self.ihi+1])) + # This is actually a pointwise norm, not quadrature'd + return numpy.sqrt(self.dx*numpy.sum(e[:, self.ilo:self.ihi+1]**2)) + +# element_norm = self.weights @ e +# return numpy.sqrt(self.dx*numpy.sum(element_norm[self.ilo:self.ihi+1]**2)) class Simulation(object): @@ -256,19 +268,19 @@ def evolve(self, num_periods=1): nx = 64 ng = 1 - g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=1) - g3 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=3) - g7 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, np=7) + g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=1) + g3 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=3) + g7 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=7) u = 1.0 - # The CFL limit for DG is reduced by a factor 1/(2 p + 1) + # The CFL limit for DG is reduced by a factor 1/(2 m + 1) s1 = Simulation(g1, u, C=0.8/(2*1+1)) s3 = Simulation(g3, u, C=0.8/(2*3+1)) s7 = Simulation(g7, u, C=0.5/(2*7+1)) # Not sure what the critical CFL is for s in [s1, s3, s7]: -# s.init_cond("sine") - s.init_cond("tophat") + s.init_cond("sine") +# s.init_cond("tophat") # Plot the initial data to show how, difference in nodal locations as # number of modes varies plot_x1, plot_a1 = g1.plotting_data() @@ -315,3 +327,32 @@ def evolve(self, num_periods=1): pyplot.ylabel(r'$a$') pyplot.legend() pyplot.show() + + # Note that the highest m (4) doesn't converge at the expected rate - + # probably limited by the time integrator. + colors = "brckgy" + symbols = "xo^ Date: Tue, 26 Feb 2019 15:45:49 +0000 Subject: [PATCH 08/33] scipy time integrator fixes all convergence issues --- advection/dg.py | 86 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 79 insertions(+), 7 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 5cd0a52..0c418d4 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -6,6 +6,7 @@ from matplotlib import pyplot import matplotlib as mpl import quadpy +from scipy.integrate import ode mpl.rcParams['mathtext.fontset'] = 'cm' mpl.rcParams['mathtext.rm'] = 'serif' @@ -257,6 +258,48 @@ def evolve(self, num_periods=1): self.t += dt + def evolve_scipy(self, num_periods=1): + """ evolve the linear advection equation using scipy """ + self.t = 0.0 + g = self.grid + + def rk_substep_scipy(t, y): + local_a = numpy.reshape(y, g.a.shape) + # Periodic BCs + local_a[:, :g.ng] = local_a[:, -2*g.ng:-g.ng] + local_a[:, -g.ng:] = local_a[:, g.ng:2*g.ng] + # Integrate flux over element + f = self.u * local_a + interior_f = g.S.T @ f + # Use Riemann solver to get fluxes between elements + al = numpy.zeros(g.nx+2*g.ng) + ar = numpy.zeros(g.nx+2*g.ng) + # i is looping over interfaces, so al is the right edge of the left + # element, etc. + for i in range(g.ilo, g.ihi+2): + al[i] = local_a[-1, i-1] + ar[i] = local_a[ 0, i ] + boundary_f = self.riemann(al, ar) + rhs = interior_f + rhs[ 0, 1:-1] += boundary_f[1:-1] + rhs[-1, 1:-1] -= boundary_f[2:] + + # Multiply by mass matrix (inverse). + rhs_i = 2 / g.dx * g.M_inv @ rhs + + return numpy.ravel(rhs_i, order='C') + + tmax = num_periods*self.period() + r = ode(rk_substep_scipy).set_integrator('dop853') + r.set_initial_value(numpy.ravel(g.a), 0) + dt = self.timestep() + + # main evolution loop + while r.successful() and r.t < tmax: + dt = min(dt, tmax - r.t) + r.integrate(r.t+dt) + g.a[:, :] = numpy.reshape(r.y, g.a.shape) + if __name__ == "__main__": @@ -265,7 +308,7 @@ def evolve(self, num_periods=1): xmin = 0.0 xmax = 1.0 - nx = 64 + nx = 8 ng = 1 g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=1) @@ -345,12 +388,41 @@ def evolve(self, num_periods=1): errs[i, j] = s.grid.norm(s.grid.a - a_init) pyplot.loglog(nxs, errs[i, :], f'{colors[i]}{symbols[i]}', label=fr'$m={{{m}}}$') - pyplot.plot(nxs, errs[0,0]/2*(nxs[0]/nxs)**2, 'b-', - label=r'$\propto (\Delta x)^2$') - pyplot.plot(nxs, errs[1,0]*2*(nxs[0]/nxs)**3, 'r-', - label=r'$\propto (\Delta x)^3$') - pyplot.plot(nxs, errs[2,0]*2*(nxs[0]/nxs)**4, 'c-', - label=r'$\propto (\Delta x)^4$') + if m < 4: + pyplot.plot(nxs, errs[i,-2]*(nxs[-2]/nxs)**(m+1), + f'{colors[i]}--', + label=fr'$\propto (\Delta x)^{{{m+1}}}$') + pyplot.xlabel(r'$N$') + pyplot.ylabel(r'$\|$Error$\|_2$') + pyplot.legend() + pyplot.show() + + + # To check that it's a limitation of the time integrator, we can use + # the scipy DOPRI8 integrator + colors = "brckgy" + symbols = "xo^ Date: Wed, 27 Feb 2019 10:03:15 +0000 Subject: [PATCH 09/33] Try using 2-norm: doesn't help? --- advection/advection.py | 11 +++++++---- advection/weno.py | 6 +++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/advection/advection.py b/advection/advection.py index 53a251f..5a69a60 100644 --- a/advection/advection.py +++ b/advection/advection.py @@ -86,14 +86,17 @@ def fill_BCs(self): # right boundary self.a[self.ihi+1+n] = self.a[self.ilo+n] - def norm(self, e): + def norm(self, e, norm="inf"): """ 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)) - return np.max(abs(e[self.ilo:self.ihi+1])) - + if norm==2: + return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2)) + elif norm=="inf": + return np.max(abs(e[self.ilo:self.ihi+1])) + else: + raise NotImplementedError("Only norms implemented are '2' and 'inf'") class Simulation(object): diff --git a/advection/weno.py b/advection/weno.py index c73550b..c0f41ff 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -405,7 +405,7 @@ def rk_substep_scipy(t, y): su.evolve(num_periods=5) - errs[-1].append(gu.norm(gu.a - ainit)) + errs[-1].append(gu.norm(gu.a - ainit, norm=2)) pyplot.clf() N = numpy.array(N, dtype=numpy.float64) @@ -468,8 +468,8 @@ def rk_substep_scipy(t, y): su.evolve_scipy(num_periods=1) # suM.evolve_scipy(num_periods=1) - errs[-1].append(gu.norm(gu.a - ainit)) -# errsM[-1].append(guM.norm(guM.a - ainit)) + errs[-1].append(gu.norm(gu.a - ainit, norm=2)) +# errsM[-1].append(guM.norm(guM.a - ainit, norm=2)) pyplot.clf() N = numpy.array(N, dtype=numpy.float64) From 6e104ba7a18e88b937972f2b42fbc3b5ba172874 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 27 Feb 2019 11:22:43 +0000 Subject: [PATCH 10/33] DG convergence plot --- advection/dg.py | 60 ++++++++++++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 23 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 0c418d4..5afc5a9 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -138,11 +138,20 @@ def init_a(x): x <= 0.666), numpy.ones_like(x), numpy.zeros_like(x)) - elif type == "sine": def init_a(x): return numpy.sin(2.0 * numpy.pi * x / (self.grid.xmax - self.grid.xmin)) + elif type == "gaussian": + def init_a(x): + local_xl = x - self.grid.dx/2 + local_xr = x + 3*self.grid.dx/2 + al = 1.0 + numpy.exp(-60.0*(local_xl - 0.5)**2) + ar = 1.0 + numpy.exp(-60.0*(local_xr - 0.5)**2) + ac = 1.0 + numpy.exp(-60.0*(x - 0.5)**2) + + return (1./6.)*(al + 4*ac + ar) + self.grid.a = init_a(self.grid.all_nodes_per_node) @@ -308,7 +317,7 @@ def rk_substep_scipy(t, y): xmin = 0.0 xmax = 1.0 - nx = 8 + nx = 4 ng = 1 g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=1) @@ -375,6 +384,7 @@ def rk_substep_scipy(t, y): # probably limited by the time integrator. colors = "brckgy" symbols = "xo^ Date: Wed, 27 Feb 2019 11:49:05 +0000 Subject: [PATCH 11/33] Fix WENO convergence plot --- advection/weno.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/advection/weno.py b/advection/weno.py index c0f41ff..cdca16d 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -379,7 +379,7 @@ def rk_substep_scipy(t, y): #-------------- RK4 - problem = "gaussian" + problem = "sine" xmin = 0.0 xmax = 1.0 @@ -400,7 +400,7 @@ def rk_substep_scipy(t, y): gu = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) su = WENOSimulation(gu, u, C=0.5, weno_order=order) - su.init_cond("gaussian") + su.init_cond("sine") ainit = su.grid.a.copy() su.evolve(num_periods=5) @@ -427,19 +427,19 @@ def rk_substep_scipy(t, y): pyplot.xlabel("N") pyplot.ylabel(r"$\| a^\mathrm{final} - a^\mathrm{init} \|_2$", fontsize=16) - pyplot.title("Convergence of Gaussian, RK4") + pyplot.title("Convergence of sine wave, RK4") pyplot.legend(frameon=False) - pyplot.savefig("weno-converge-gaussian-rk4.pdf") + pyplot.savefig("weno-converge-sine-rk4.pdf") # pyplot.show() -#-------------- Gaussian +#-------------- Sine wave, 8th order time integrator - problem = "gaussian" + problem = "sine" xmin = 0.0 xmax = 1.0 - orders = [3, 4, 5, 6] + orders = [3, 5, 7] N = [24, 32, 54, 64, 81, 108, 128] # N = [32, 64, 108, 128] @@ -461,8 +461,8 @@ def rk_substep_scipy(t, y): # guM = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) # suM = WENOMSimulation(guM, u, C=0.5, weno_order=order) - su.init_cond("gaussian") -# suM.init_cond("gaussian") + su.init_cond("sine") +# suM.init_cond("sine") ainit = su.grid.a.copy() su.evolve_scipy(num_periods=1) @@ -480,9 +480,9 @@ def rk_substep_scipy(t, y): # pyplot.scatter(N, errsM[n_order], # color=colors[n_order], # label=r"WENOM, $r={}$".format(order)) - pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-2), + pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-1), linestyle="--", color=colors[n_order], - label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-2)) + label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-1)) # pyplot.plot(N, errs[n_order][len(N)-1]*(N[len(N)-1]/N)**4, # color="k", label=r"$\mathcal{O}(\Delta x^4)$") @@ -494,9 +494,10 @@ def rk_substep_scipy(t, y): pyplot.xlabel("N") pyplot.ylabel(r"$\| a^\mathrm{final} - a^\mathrm{init} \|_2$", fontsize=16) - pyplot.title("Convergence of Gaussian, DOPRK8") + pyplot.title("Convergence of sine wave, DOPRK8") - pyplot.legend(frameon=False) - pyplot.savefig("weno-converge-gaussian.pdf") + lgd = ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) + pyplot.savefig("weno-converge-sine.pdf", + bbox_extra_artists=(lgd,), bbox_inches='tight') # pyplot.show() \ No newline at end of file From afb7f9f5049a2ddbb31013e4eae9c7989a2bcaab Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 27 Feb 2019 11:54:56 +0000 Subject: [PATCH 12/33] Switch RK4 case back to Gaussian --- advection/weno.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/advection/weno.py b/advection/weno.py index cdca16d..4f66faf 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -379,7 +379,7 @@ def rk_substep_scipy(t, y): #-------------- RK4 - problem = "sine" + problem = "gaussian" xmin = 0.0 xmax = 1.0 @@ -400,7 +400,7 @@ def rk_substep_scipy(t, y): gu = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) su = WENOSimulation(gu, u, C=0.5, weno_order=order) - su.init_cond("sine") + su.init_cond("gaussian") ainit = su.grid.a.copy() su.evolve(num_periods=5) @@ -427,10 +427,10 @@ def rk_substep_scipy(t, y): pyplot.xlabel("N") pyplot.ylabel(r"$\| a^\mathrm{final} - a^\mathrm{init} \|_2$", fontsize=16) - pyplot.title("Convergence of sine wave, RK4") + pyplot.title("Convergence of Gaussian, RK4") pyplot.legend(frameon=False) - pyplot.savefig("weno-converge-sine-rk4.pdf") + pyplot.savefig("weno-converge-gaussian-rk4.pdf") # pyplot.show() #-------------- Sine wave, 8th order time integrator From eea5a37d3bc6a66442d20599c109401aee7f4aa6 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 27 Feb 2019 16:45:50 +0000 Subject: [PATCH 13/33] Show modal structure within cells --- advection/dg.py | 84 ++++++++++++++++++------------------------------- 1 file changed, 31 insertions(+), 53 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 5afc5a9..d8d2531 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -34,7 +34,7 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): self.dx = (xmax - xmin)/(nx) self.x = xmin + (numpy.arange(nx+2*ng)-ng+0.5)*self.dx self.xl = xmin + (numpy.arange(nx+2*ng)-ng)*self.dx - self.xr = xmin + (numpy.arange(nx+2*ng)+1.0)*self.dx + self.xr = xmin + (numpy.arange(nx+2*ng)-ng+1.0)*self.dx # storage for the solution # These are the modes of the solution at each point, so the @@ -87,6 +87,21 @@ def nodal_to_modal(self): def plotting_data(self): return (self.all_nodes, self.a.ravel(order='F')) + + def plotting_data_high_order(self, npoints=50): + assert npoints > 2 + p_nodes = numpy.zeros(npoints*(self.nx+2*self.ng), dtype=numpy.float64) + p_data = numpy.zeros_like(p_nodes) + for i in range(self.nx+2*self.ng): + p_nodes[npoints*i:npoints*(i+1)] = numpy.linspace(self.xl[i], + self.xr[i], + npoints) + modal = self.V_inv @ self.a[:, i] + for p in range(self.m+1): + modal[p] /= numpy.sqrt(2/(2*p+1)) + scaled_x = 2 * (p_nodes[npoints*i:npoints*(i+1)] - self.x[i]) / self.dx + p_data[npoints*i:npoints*(i+1)] = numpy.polynomial.legendre.legval(scaled_x, modal) + return p_nodes, p_data def scratch_array(self): """ return a scratch array dimensioned for our grid """ @@ -145,7 +160,7 @@ def init_a(x): elif type == "gaussian": def init_a(x): local_xl = x - self.grid.dx/2 - local_xr = x + 3*self.grid.dx/2 + local_xr = x + self.grid.dx/2 al = 1.0 + numpy.exp(-60.0*(local_xl - 0.5)**2) ar = 1.0 + numpy.exp(-60.0*(local_xr - 0.5)**2) ac = 1.0 + numpy.exp(-60.0*(x - 0.5)**2) @@ -313,72 +328,35 @@ def rk_substep_scipy(t, y): if __name__ == "__main__": # ------------------------------------------------------------------------- - # DG of sine wave + # Show the "grid" using a sine wave xmin = 0.0 xmax = 1.0 nx = 4 ng = 1 - g1 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=1) - g3 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=3) - g7 = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=7) - u = 1.0 - # The CFL limit for DG is reduced by a factor 1/(2 m + 1) - s1 = Simulation(g1, u, C=0.8/(2*1+1)) - s3 = Simulation(g3, u, C=0.8/(2*3+1)) - s7 = Simulation(g7, u, C=0.5/(2*7+1)) # Not sure what the critical CFL is - for s in [s1, s3, s7]: + colors="kbr" + symbols="sox" + for m in range(1, 4): + g = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=m) + s = Simulation(g, u, C=0.5/(2*m+1)) s.init_cond("sine") -# s.init_cond("tophat") - # Plot the initial data to show how, difference in nodal locations as - # number of modes varies - plot_x1, plot_a1 = g1.plotting_data() - a1init = plot_a1.copy() - plot_x3, plot_a3 = g3.plotting_data() - a3init = plot_a3.copy() - plot_x7, plot_a7 = g7.plotting_data() - a7init = plot_a7.copy() - pyplot.plot(plot_x1, plot_a1, 'k>', label=r"$p=1$") - pyplot.plot(plot_x3, plot_a3, 'bo', label=r"$p=3$") - pyplot.plot(plot_x7, plot_a7, 'r^', label=r"$p=7$") - pyplot.xlim(0, 1) - pyplot.xlabel(r'$x$') - pyplot.ylabel(r'$a$') - pyplot.legend() - pyplot.show() - - s1.evolve(num_periods=1) - plot_x1, plot_a1 = g1.plotting_data() - pyplot.plot(plot_x1, plot_a1, 'k>', label=r"$p=1$") - pyplot.plot(plot_x1, a1init, 'r^', label=r"$p=1$, initial") - pyplot.xlim(0, 1) - pyplot.xlabel(r'$x$') - pyplot.ylabel(r'$a$') - pyplot.legend() - pyplot.show() - - s3.evolve(num_periods=1) - plot_x3, plot_a3 = g3.plotting_data() - pyplot.plot(plot_x3, plot_a3, 'bo', label=r"$p=3$") - pyplot.plot(plot_x3, a3init, 'r^', label=r"$p=3$, initial") + plot_x, plot_a = g.plotting_data() + plot_x_hi, plot_a_hi = g.plotting_data_high_order() + pyplot.plot(plot_x, plot_a, f'{colors[m-1]}{symbols[m-1]}', + label=fr"Nodes, $m={{{m}}}$") + pyplot.plot(plot_x_hi, plot_a_hi, f'{colors[m-1]}:', + label=fr"Modes, $m={{{m}}}$") pyplot.xlim(0, 1) + pyplot.vlines([0.25, 0.5, 0.75], -2, 2, linestyles='--') + pyplot.ylim(-1, 1) pyplot.xlabel(r'$x$') pyplot.ylabel(r'$a$') pyplot.legend() pyplot.show() - s7.evolve(num_periods=1) - plot_x7, plot_a7 = g7.plotting_data() - pyplot.plot(plot_x7, plot_a7, 'bo', label=r"$p=7$") - pyplot.plot(plot_x7, a7init, 'r^', label=r"$p=7$, initial") - pyplot.xlim(0, 1) - pyplot.xlabel(r'$x$') - pyplot.ylabel(r'$a$') - pyplot.legend() - pyplot.show() # Note that the highest m (4) doesn't converge at the expected rate - # probably limited by the time integrator. From 1ee4380b9a23ffaaa95fdc1f259f69846eaf0bd9 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Fri, 1 Mar 2019 14:54:44 +0000 Subject: [PATCH 14/33] Moment limiting. Not correct yet --- advection/dg.py | 156 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 129 insertions(+), 27 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index d8d2531..03fdea9 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -72,11 +72,10 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): self.all_nodes_per_node[:, i] = (self.x[i] + self.nodes * self.dx / 2) -# def modal_to_nodal(self): -# nodal = numpy.zeros_like(self.a) -# for i in range(self.nx+2*self.ng): -# nodal[:, i] = self.V @ self.a[:, i] -# return nodal + def modal_to_nodal(self, modal): + for i in range(self.nx+2*self.ng): + self.a[:, i] = self.V @ modal[:, i] + return self.a def nodal_to_modal(self): modal = numpy.zeros_like(self.a) @@ -137,13 +136,30 @@ def norm(self, e): # return numpy.sqrt(self.dx*numpy.sum(element_norm[self.ilo:self.ihi+1]**2)) +def minmod3(a1, a2, a3): + """ + Utility function that does minmod on three inputs + """ + signs1 = a1 * a2 + signs2 = a1 * a3 + signs3 = a2 * a3 + same_sign = numpy.logical_and(numpy.logical_and(signs1 > 0, + signs2 > 0), + signs3 > 0) + minmod = numpy.min(numpy.abs(numpy.vstack((a1, a2, a3))), + axis=0) * numpy.sign(a1) + + return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) + + class Simulation(object): - def __init__(self, grid, u, C=0.8): + def __init__(self, grid, u, C=0.8, limiter=None): self.grid = grid self.t = 0.0 # simulation time self.u = u # the constant advective velocity self.C = C # CFL number + self.limiter = limiter # What it says. def init_cond(self, type="sine"): """ initialize the data """ @@ -184,6 +200,8 @@ def states(self): # Evaluate the nodal values at the domain edges g = self.grid + # Extract nodal values + al = numpy.zeros(g.nx+2*g.ng) ar = numpy.zeros(g.nx+2*g.ng) @@ -193,30 +211,70 @@ def states(self): al[i] = g.a[-1, i-1] ar[i] = g.a[ 0, i ] - # TODO: remove this return, fix the limiting. return al, ar - # Limiting, using minmod (Hesthaven p 443-4) - theta = 2 - a_modal = g.nodal_to_modal() - a_average = a_modal[0, :] - a_m = numpy.zeros_like(a_average) - a_m[1:] = a_average[:-1] - a_el = g.a[0 , :] - a_er = g.a[-1, :] - - check_left = a_average - minmod([a_average - a_el, - a_average - a_m, - a_er - a_average]) - check_right = a_average + minmod([a_average - a_el, - a_average - m, - a_er - a_average]) - ids = numpy.logical_or(numpy.isclose(a_el, check_left), - numpy.isclose(a_er, check_right)) - + def limit(self): + """ + After evolution, limit the slopes. + """ - return al, ar + # Evaluate the nodal values at the domain edges + g = self.grid + # Limiting! + + if self.limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[1:, :] = 0 + a_cell = (g.V @ a_zeromode)[0,:] + a_minus = g.a[0, :] + a_plus = g.a[-1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[i+1, 1:-1] + a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) + a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor +# print("Original", i, a_modal[i+1, :]) +# print("Updated", i, updated_mode) + did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], + updated_mode[1:-1]) + a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] +# print("Limited", i, a_modal[i+1, :]) + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + did_it_limit) + # Get back nodal values + g.a = g.modal_to_nodal(a_modal) + + return None + def riemann(self, al, ar): """ Riemann problem for advection -- this is simply upwinding, @@ -273,12 +331,15 @@ def evolve(self, num_periods=1): a_start = g.a.copy() k1 = dt * self.rk_substep() g.a = a_start + k1 + self.limit() a1 = g.a.copy() k2 = dt * self.rk_substep() g.a = (3 * a_start + a1 + k2) / 4 + self.limit() a2 = g.a.copy() k3 = dt * self.rk_substep() g.a = (a_start + 2 * a2 + 2 * k3) / 3 + self.limit() self.t += dt @@ -322,11 +383,52 @@ def rk_substep_scipy(t, y): while r.successful() and r.t < tmax: dt = min(dt, tmax - r.t) r.integrate(r.t+dt) - g.a[:, :] = numpy.reshape(r.y, g.a.shape) + g.a[:, :] = numpy.reshape(r.y, g.a.shape) + self.limit() + r.y = numpy.ravel(g.a) if __name__ == "__main__": + # Checking the limiter + g = Grid1d(16, 1, 0, 1, 3) + s = Simulation(g, 1, 0.5/7, "moment") + s.init_cond("sine") + s.limit() + + + # Runs with limiter + g_sin_nolimit = Grid1d(16, 1, 0, 1, 3) + g_hat_nolimit = Grid1d(16, 1, 0, 1, 3) + g_sin_moment = Grid1d(16, 1, 0, 1, 3) + g_hat_moment = Grid1d(16, 1, 0, 1, 3) + s_sin_nolimit = Simulation(g_sin_nolimit, 1, 0.5/7, limiter=None) + s_hat_nolimit = Simulation(g_hat_nolimit, 1, 0.5/7, limiter=None) + s_sin_moment = Simulation(g_sin_moment, 1, 0.5/7, limiter="moment") + s_hat_moment = Simulation(g_hat_moment, 1, 0.5/7, limiter="moment") + for s in s_sin_nolimit, s_sin_moment: + s.init_cond("sine") + for s in s_hat_nolimit, s_hat_moment: + s.init_cond("tophat") + for s in s_sin_nolimit, s_sin_moment, s_hat_nolimit, s_hat_moment: + s.evolve() + fig, axes = pyplot.subplots(2, 2) + for i, g in enumerate((g_sin_nolimit, g_hat_nolimit)): + plot_x, plot_a = g.plotting_data() + axes[0, i].plot(plot_x, plot_a) + axes[0, i].set_xlim(0, 1) + axes[0, i].set_xlabel(r"$x$") + axes[0, i].set_ylabel(r"$a$") + for i, g in enumerate((g_sin_moment, g_hat_moment)): + plot_x, plot_a = g.plotting_data() + axes[1, i].plot(plot_x, plot_a) + axes[1, i].set_xlim(0, 1) + axes[1, i].set_xlabel(r"$x$") + axes[1, i].set_ylabel(r"$a$") + fig.tight_layout() + pyplot.show() + 1/0 + # ------------------------------------------------------------------------- # Show the "grid" using a sine wave From 0f054d6181a0aa97dd02c51eb7e14276a11abc40 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Fri, 1 Mar 2019 17:04:47 +0000 Subject: [PATCH 15/33] Limiting works, but broke scipy evolution --- advection/dg.py | 81 ++++++++++++++++++++++++++++++++++--------------- 1 file changed, 57 insertions(+), 24 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 03fdea9..8aa104b 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -269,11 +269,11 @@ def limit(self): a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] # print("Limited", i, a_modal[i+1, :]) limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], - did_it_limit) + numpy.logical_not(did_it_limit)) # Get back nodal values g.a = g.modal_to_nodal(a_modal) - return None + return g.a def riemann(self, al, ar): """ @@ -331,15 +331,15 @@ def evolve(self, num_periods=1): a_start = g.a.copy() k1 = dt * self.rk_substep() g.a = a_start + k1 - self.limit() + g.a = self.limit() a1 = g.a.copy() k2 = dt * self.rk_substep() g.a = (3 * a_start + a1 + k2) / 4 - self.limit() + g.a = self.limit() a2 = g.a.copy() k3 = dt * self.rk_substep() g.a = (a_start + 2 * a2 + 2 * k3) / 3 - self.limit() + g.a = self.limit() self.t += dt @@ -390,13 +390,6 @@ def rk_substep_scipy(t, y): if __name__ == "__main__": - # Checking the limiter - g = Grid1d(16, 1, 0, 1, 3) - s = Simulation(g, 1, 0.5/7, "moment") - s.init_cond("sine") - s.limit() - - # Runs with limiter g_sin_nolimit = Grid1d(16, 1, 0, 1, 3) g_hat_nolimit = Grid1d(16, 1, 0, 1, 3) @@ -410,24 +403,28 @@ def rk_substep_scipy(t, y): s.init_cond("sine") for s in s_hat_nolimit, s_hat_moment: s.init_cond("tophat") + x_sin_start, a_sin_start = g_sin_nolimit.plotting_data() + x_hat_start, a_hat_start = g_hat_nolimit.plotting_data() for s in s_sin_nolimit, s_sin_moment, s_hat_nolimit, s_hat_moment: s.evolve() - fig, axes = pyplot.subplots(2, 2) + fig, axes = pyplot.subplots(1, 2) + axes[0].plot(x_sin_start, a_sin_start, 'k-', label='Exact') + axes[1].plot(x_hat_start, a_hat_start, 'k-', label='Exact') for i, g in enumerate((g_sin_nolimit, g_hat_nolimit)): - plot_x, plot_a = g.plotting_data() - axes[0, i].plot(plot_x, plot_a) - axes[0, i].set_xlim(0, 1) - axes[0, i].set_xlabel(r"$x$") - axes[0, i].set_ylabel(r"$a$") + axes[i].plot(*g.plotting_data(), 'b--', label='No limiting') + axes[i].set_xlim(0, 1) + axes[i].set_xlabel(r"$x$") + axes[i].set_ylabel(r"$a$") for i, g in enumerate((g_sin_moment, g_hat_moment)): - plot_x, plot_a = g.plotting_data() - axes[1, i].plot(plot_x, plot_a) - axes[1, i].set_xlim(0, 1) - axes[1, i].set_xlabel(r"$x$") - axes[1, i].set_ylabel(r"$a$") + axes[i].plot(*g.plotting_data(), 'r:', label='Moment limiting') + axes[i].set_xlim(0, 1) + axes[i].set_xlabel(r"$x$") + axes[i].set_ylabel(r"$a$") fig.tight_layout() + lgd = axes[1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) +# fig.savefig('dg_limiter.pdf', +# bbox_extra_artists=(lgd,), bbox_inches='tight') pyplot.show() - 1/0 # ------------------------------------------------------------------------- # Show the "grid" using a sine wave @@ -521,4 +518,40 @@ def rk_substep_scipy(t, y): fig.savefig('dg_convergence_sine.pdf', bbox_extra_artists=(lgd,), bbox_inches='tight') pyplot.show() + + + + # Check convergence with the limiter on + colors = "brckgy" + symbols = "xo^ Date: Sat, 2 Mar 2019 12:36:37 +0000 Subject: [PATCH 16/33] Fix various PEP8 issues --- advection/dg.py | 114 +++++++++++++++++++++--------------------------- 1 file changed, 50 insertions(+), 64 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 8aa104b..14f1fcd 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -3,6 +3,7 @@ """ import numpy +from numpy.polynomial import legendre from matplotlib import pyplot import matplotlib as mpl import quadpy @@ -47,7 +48,7 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): self.nodes = GL.points self.weights = GL.weights # To go from modal to nodal we need the Vandermonde matrix - self.V = numpy.polynomial.legendre.legvander(self.nodes, m) + self.V = legendre.legvander(self.nodes, m) c = numpy.eye(m+1) # Orthonormalize for p in range(m+1): @@ -57,8 +58,8 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): self.M = numpy.linalg.inv(self.V @ self.V.T) self.M_inv = self.V @ self.V.T # Derivatives of Legendre polynomials lead to derivatives of V - dV = numpy.polynomial.legendre.legval(self.nodes, - numpy.polynomial.legendre.legder(c)).T + dV = legendre.legval(self.nodes, + legendre.legder(c)).T self.D = dV @ self.V_inv # Stiffness matrix for the interior flux self.S = self.M @ self.D @@ -68,7 +69,7 @@ def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): self.all_nodes_per_node = numpy.zeros_like(self.a) for i in range(nx+2*ng): self.all_nodes[(m+1)*i:(m+1)*(i+1)] = (self.x[i] + - self.nodes * self.dx / 2) + self.nodes * self.dx / 2) self.all_nodes_per_node[:, i] = (self.x[i] + self.nodes * self.dx / 2) @@ -86,7 +87,7 @@ def nodal_to_modal(self): def plotting_data(self): return (self.all_nodes, self.a.ravel(order='F')) - + def plotting_data_high_order(self, npoints=50): assert npoints > 2 p_nodes = numpy.zeros(npoints*(self.nx+2*self.ng), dtype=numpy.float64) @@ -98,8 +99,9 @@ def plotting_data_high_order(self, npoints=50): modal = self.V_inv @ self.a[:, i] for p in range(self.m+1): modal[p] /= numpy.sqrt(2/(2*p+1)) - scaled_x = 2 * (p_nodes[npoints*i:npoints*(i+1)] - self.x[i]) / self.dx - p_data[npoints*i:npoints*(i+1)] = numpy.polynomial.legendre.legval(scaled_x, modal) + scaled_x = 2 * (p_nodes[npoints*i:npoints*(i+1)] - + self.x[i]) / self.dx + p_data[npoints*i:npoints*(i+1)] = legendre.legval(scaled_x, modal) return p_nodes, p_data def scratch_array(self): @@ -116,9 +118,9 @@ def fill_BCs(self): self.a[:, self.ihi+1+n] = self.a[:, self.ilo+n] def norm(self, e): - """ + """ Return the norm of quantity e which lives on the grid. - + This is the 'broken norm': the quantity is integrated over each individual element using Gauss-Lobatto quadrature (as we have those nodes and weights), and the 2-norm of the result is then returned. @@ -126,15 +128,9 @@ def norm(self, e): if not numpy.allclose(e.shape, self.all_nodes_per_node.shape): return None - - # This is L_inf norm... -# return numpy.max(abs(e[:, self.ilo:self.ihi+1])) # This is actually a pointwise norm, not quadrature'd return numpy.sqrt(self.dx*numpy.sum(e[:, self.ilo:self.ihi+1]**2)) -# element_norm = self.weights @ e -# return numpy.sqrt(self.dx*numpy.sum(element_norm[self.ilo:self.ihi+1]**2)) - def minmod3(a1, a2, a3): """ @@ -148,7 +144,7 @@ def minmod3(a1, a2, a3): signs3 > 0) minmod = numpy.min(numpy.abs(numpy.vstack((a1, a2, a3))), axis=0) * numpy.sign(a1) - + return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) @@ -180,9 +176,8 @@ def init_a(x): al = 1.0 + numpy.exp(-60.0*(local_xl - 0.5)**2) ar = 1.0 + numpy.exp(-60.0*(local_xr - 0.5)**2) ac = 1.0 + numpy.exp(-60.0*(x - 0.5)**2) - - return (1./6.)*(al + 4*ac + ar) + return (1./6.)*(al + 4*ac + ar) self.grid.a = init_a(self.grid.all_nodes_per_node) @@ -201,7 +196,7 @@ def states(self): g = self.grid # Extract nodal values - + al = numpy.zeros(g.nx+2*g.ng) ar = numpy.zeros(g.nx+2*g.ng) @@ -209,10 +204,10 @@ def states(self): # element, etc. for i in range(g.ilo, g.ihi+2): al[i] = g.a[-1, i-1] - ar[i] = g.a[ 0, i ] + ar[i] = g.a[0, i] return al, ar - + def limit(self): """ After evolution, limit the slopes. @@ -222,9 +217,9 @@ def limit(self): g = self.grid # Limiting! - + if self.limiter == "moment": - + # Limiting, using moment limiting (Hesthaven p 445-7) theta = 2 a_modal = g.nodal_to_modal() @@ -235,19 +230,19 @@ def limit(self): # Get the cell average and the nodal values at the boundaries a_zeromode = a_modal.copy() a_zeromode[1:, :] = 0 - a_cell = (g.V @ a_zeromode)[0,:] + a_cell = (g.V @ a_zeromode)[0, :] a_minus = g.a[0, :] a_plus = g.a[-1, :] - # From the cell averages and boundary values we can construct + # From the cell averages and boundary values we can construct # alternate values at the boundaries a_left = numpy.zeros(g.nx+2*g.ng) a_right = numpy.zeros(g.nx+2*g.ng) a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], - a_cell[1:-1] - a_cell[:-2], - a_cell[2:] - a_cell[1:-1]) - a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], a_cell[1:-1] - a_cell[:-2], a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) limiting_todo[1:-1] = numpy.logical_not( numpy.logical_and(numpy.isclose(a_left[1:-1], a_minus[1:-1]), @@ -262,19 +257,16 @@ def limit(self): a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) updated_mode[1:-1] = minmod3(a1, a2, a3) / factor -# print("Original", i, a_modal[i+1, :]) -# print("Updated", i, updated_mode) did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], updated_mode[1:-1]) a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] -# print("Limited", i, a_modal[i+1, :]) limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], - numpy.logical_not(did_it_limit)) + numpy.logical_not(did_it_limit)) # Get back nodal values g.a = g.modal_to_nodal(a_modal) - + return g.a - + def riemann(self, al, ar): """ Riemann problem for advection -- this is simply upwinding, @@ -300,7 +292,7 @@ def rk_substep(self): # Use Riemann solver to get fluxes between elements boundary_f = self.riemann(*self.states()) rhs = interior_f - rhs[ 0, 1:-1] += boundary_f[1:-1] + rhs[0, 1:-1] += boundary_f[1:-1] rhs[-1, 1:-1] -= boundary_f[2:] # Multiply by mass matrix (inverse). @@ -347,7 +339,7 @@ def evolve_scipy(self, num_periods=1): """ evolve the linear advection equation using scipy """ self.t = 0.0 g = self.grid - + def rk_substep_scipy(t, y): local_a = numpy.reshape(y, g.a.shape) # Periodic BCs @@ -363,15 +355,15 @@ def rk_substep_scipy(t, y): # element, etc. for i in range(g.ilo, g.ihi+2): al[i] = local_a[-1, i-1] - ar[i] = local_a[ 0, i ] + ar[i] = local_a[0, i] boundary_f = self.riemann(al, ar) rhs = interior_f - rhs[ 0, 1:-1] += boundary_f[1:-1] + rhs[0, 1:-1] += boundary_f[1:-1] rhs[-1, 1:-1] -= boundary_f[2:] - + # Multiply by mass matrix (inverse). rhs_i = 2 / g.dx * g.M_inv @ rhs - + return numpy.ravel(rhs_i, order='C') tmax = num_periods*self.period() @@ -383,10 +375,9 @@ def rk_substep_scipy(t, y): while r.successful() and r.t < tmax: dt = min(dt, tmax - r.t) r.integrate(r.t+dt) - g.a[:, :] = numpy.reshape(r.y, g.a.shape) - self.limit() - r.y = numpy.ravel(g.a) - + g.a[:, :] = numpy.reshape(r.y, g.a.shape) + self.limit() + if __name__ == "__main__": @@ -422,10 +413,10 @@ def rk_substep_scipy(t, y): axes[i].set_ylabel(r"$a$") fig.tight_layout() lgd = axes[1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) -# fig.savefig('dg_limiter.pdf', +# fig.savefig('dg_limiter.pdf', # bbox_extra_artists=(lgd,), bbox_inches='tight') pyplot.show() - + # ------------------------------------------------------------------------- # Show the "grid" using a sine wave @@ -436,8 +427,8 @@ def rk_substep_scipy(t, y): u = 1.0 - colors="kbr" - symbols="sox" + colors = "kbr" + symbols = "sox" for m in range(1, 4): g = Grid1d(nx, ng, xmin=xmin, xmax=xmax, m=m) s = Simulation(g, u, C=0.5/(2*m+1)) @@ -456,8 +447,7 @@ def rk_substep_scipy(t, y): pyplot.legend() pyplot.show() - - # Note that the highest m (4) doesn't converge at the expected rate - + # Note that the highest m (4) doesn't converge at the expected rate - # probably limited by the time integrator. colors = "brckgy" symbols = "xo^ Date: Sat, 2 Mar 2019 17:11:25 +0000 Subject: [PATCH 17/33] Limiting works with scipy ode --- advection/dg.py | 47 ++++++++++++++++++++++++++++------------------- 1 file changed, 28 insertions(+), 19 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 14f1fcd..e6f1bc8 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -367,16 +367,24 @@ def rk_substep_scipy(t, y): return numpy.ravel(rhs_i, order='C') tmax = num_periods*self.period() - r = ode(rk_substep_scipy).set_integrator('dop853') - r.set_initial_value(numpy.ravel(g.a), 0) - dt = self.timestep() # main evolution loop - while r.successful() and r.t < tmax: - dt = min(dt, tmax - r.t) + while self.t < tmax: + # fill the boundary conditions + g.fill_BCs() + + # get the timestep + dt = self.timestep() + + if self.t + dt > tmax: + dt = tmax - self.t + + r = ode(rk_substep_scipy).set_integrator('dop853') + r.set_initial_value(numpy.ravel(g.a), self.t) r.integrate(r.t+dt) - g.a[:, :] = numpy.reshape(r.y, g.a.shape) - self.limit() + g.a[:, :] = numpy.reshape(r.y, g.a.shape) + g.a = self.limit() + self.t += dt if __name__ == "__main__": @@ -511,6 +519,7 @@ def rk_substep_scipy(t, y): # Check convergence with the limiter on colors = "brckgy" symbols = "xo^ Date: Sun, 3 Mar 2019 17:41:13 +0000 Subject: [PATCH 18/33] Better convergence figures --- advection/dg.py | 56 ++++++++++++++++++++----------------------------- 1 file changed, 23 insertions(+), 33 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index e6f1bc8..33fe81e 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -459,7 +459,7 @@ def rk_substep_scipy(t, y): # probably limited by the time integrator. colors = "brckgy" symbols = "xo^ Date: Mon, 4 Mar 2019 13:18:19 +0000 Subject: [PATCH 19/33] Figure of runtime efficiency --- advection/compare_schemes.py | 115 +++++++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 advection/compare_schemes.py diff --git a/advection/compare_schemes.py b/advection/compare_schemes.py new file mode 100644 index 0000000..f1c379b --- /dev/null +++ b/advection/compare_schemes.py @@ -0,0 +1,115 @@ +""" +Comparing efficiency (accuracy vs runtime and memory) for various schemes. + +Note: needs 'pip install memory_profiler' first. + +At least for now, the memory usage seems to be completely misleading. +""" + +import timeit +from memory_profiler import memory_usage +import numpy +from matplotlib import pyplot +import advection +import weno +import dg + +def run_weno(order=2, nx=32): + g = advection.Grid1d(nx, ng=order+1, xmin=0, xmax=1) + s = weno.WENOSimulation(g, 1, C=0.5, weno_order=order) + s.init_cond("sine") + s.evolve_scipy() + +def run_dg(m=3, nx=8, limiter=None): + g = dg.Grid1d(nx, ng=1, xmin=0, xmax=1, m=m) + s = dg.Simulation(g, 1, C=0.5/(2*m+1), limiter=limiter) + s.init_cond("sine") + s.evolve_scipy() + +def weno_string(order=2, nx=32): + return f"run_weno(order={order}, nx={nx})" + +def dg_string(m=3, nx=8, limiter=None): + return f"run_dg(m={m}, nx={nx}, limiter='{limiter}')" + +def time_weno(order=2, nx=32): + return timeit.timeit(weno_string(order, nx), + globals=globals(), number=1) + +def time_dg(m=3, nx=8, limiter=None): + return timeit.timeit(dg_string(m, nx, limiter), + globals=globals(), number=1) + +def errs_weno(order=2, nx=32): + g = advection.Grid1d(nx, ng=order+1, xmin=0, xmax=1) + s = weno.WENOSimulation(g, 1, C=0.5, weno_order=order) + s.init_cond("sine") + a_init = g.a.copy() + s.evolve_scipy() + return g.norm(g.a - a_init) + +def errs_dg(m=3, nx=8, limiter=None): + g = dg.Grid1d(nx, ng=1, xmin=0, xmax=1, m=m) + s = dg.Simulation(g, 1, C=0.5/(2*m+1), limiter=limiter) + s.init_cond("sine") + a_init = g.a.copy() + s.evolve_scipy() + return g.norm(g.a - a_init) + +weno_orders = [3, 5, 7] +weno_N = [8, 12, 16, 24, 32, 54, 64] +#weno_orders = [3] +#weno_N = [24, 32, 54] +weno_times = numpy.zeros((len(weno_orders), len(weno_N))) +weno_errs = numpy.zeros_like(weno_times) +weno_memory = numpy.zeros_like(weno_times) + +for i_o, order in enumerate(weno_orders): + for i_n, nx in enumerate(weno_N): + print(f"WENO{order}, {nx} points") + weno_errs[i_o, i_n] = errs_weno(order, nx) + weno_times[i_o, i_n] = time_weno(order, nx) + weno_memory[i_o, i_n] = max(memory_usage((run_weno, (order, nx)))) + +dg_ms = [2, 3, 4] +dg_N = 2**numpy.array(range(3, 9)) +#dg_ms = [2] +#dg_N = 2**numpy.array(range(3, 6)) +dg_nolimit_times = numpy.zeros((len(dg_ms), len(dg_N))) +dg_moment_times = numpy.zeros((len(dg_ms), len(dg_N))) +dg_nolimit_errs = numpy.zeros_like(dg_nolimit_times) +dg_moment_errs = numpy.zeros_like(dg_moment_times) +dg_nolimit_memory = numpy.zeros_like(dg_nolimit_times) +dg_moment_memory = numpy.zeros_like(dg_moment_times) + +for i_m, m in enumerate(dg_ms): + for i_n, nx in enumerate(dg_N): + print(f"DG, m={m}, {nx} points") + dg_nolimit_errs[i_m, i_n] = errs_dg(m, nx, None) + dg_nolimit_times[i_m, i_n] = time_dg(m, nx, None) + dg_nolimit_memory[i_m, i_n] = max(memory_usage((run_dg, + (m, nx, None)))) + dg_moment_errs[i_m, i_n] = errs_dg(m, nx, 'moment') + dg_moment_times[i_m, i_n] = time_dg(m, nx, 'moment') + dg_moment_memory[i_m, i_n] = max(memory_usage((run_dg, + (m, nx, 'moment')))) + +colors = "brk" +fig, ax = pyplot.subplots(1, 1) +for i_o, order in enumerate(weno_orders): + ax.loglog(weno_times[i_o, :], weno_errs[i_o, :], f'{colors[i_o]}o-', + label=f'WENO, order={order}') +colors = "gyc" +for i_m, m in enumerate(dg_ms): + ax.loglog(dg_nolimit_times[i_m, :], dg_nolimit_errs[i_m, :], + f'{colors[i_m]}s--', label=rf'DG, no limiter, $m={m}$') + ax.loglog(dg_moment_times[i_m, :], dg_moment_errs[i_m, :], + f'{colors[i_m]}^:', label=rf'DG, moment limiter, $m={m}$') +ax.set_xlabel("Runtime [s]") +ax.set_ylabel(r"$\|$Error$\|_2$") +ax.set_title("Efficiency of WENO vs DG") +fig.tight_layout() +lgd = ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) +fig.savefig('dg_weno_efficiency.pdf', + bbox_extra_artists=(lgd,), bbox_inches='tight') +pyplot.show() From 1f93f402cd9fff1c739a77b4e8d0c7b186d5aad1 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Mon, 4 Mar 2019 13:22:27 +0000 Subject: [PATCH 20/33] Output earlier DG figures --- advection/dg.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index 33fe81e..327c858 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -421,8 +421,8 @@ def rk_substep_scipy(t, y): axes[i].set_ylabel(r"$a$") fig.tight_layout() lgd = axes[1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) -# fig.savefig('dg_limiter.pdf', -# bbox_extra_artists=(lgd,), bbox_inches='tight') + fig.savefig('dg_limiter.pdf', + bbox_extra_artists=(lgd,), bbox_inches='tight') pyplot.show() # ------------------------------------------------------------------------- @@ -453,6 +453,7 @@ def rk_substep_scipy(t, y): pyplot.xlabel(r'$x$') pyplot.ylabel(r'$a$') pyplot.legend() + fig.savefig("dg_grid.pdf") pyplot.show() # Note that the highest m (4) doesn't converge at the expected rate - From b3af182c9db1bd3e165f1ef15d0b197562fc7cac Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Mon, 4 Mar 2019 13:29:30 +0000 Subject: [PATCH 21/33] Fix figure output --- advection/dg.py | 1 + 1 file changed, 1 insertion(+) diff --git a/advection/dg.py b/advection/dg.py index 327c858..cc33370 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -435,6 +435,7 @@ def rk_substep_scipy(t, y): u = 1.0 + pyplot.clf() colors = "kbr" symbols = "sox" for m in range(1, 4): From cf000b91c4f0496246963b6f4c2c96b94952b38e Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Mon, 4 Mar 2019 16:54:08 +0000 Subject: [PATCH 22/33] Fix figure output --- advection/dg.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/advection/dg.py b/advection/dg.py index cc33370..6b02ca7 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -435,7 +435,6 @@ def rk_substep_scipy(t, y): u = 1.0 - pyplot.clf() colors = "kbr" symbols = "sox" for m in range(1, 4): @@ -454,7 +453,7 @@ def rk_substep_scipy(t, y): pyplot.xlabel(r'$x$') pyplot.ylabel(r'$a$') pyplot.legend() - fig.savefig("dg_grid.pdf") + pyplot.savefig("dg_grid.pdf") pyplot.show() # Note that the highest m (4) doesn't converge at the expected rate - From a480adbc3379ff9238f9ab539b6e815672d1bb65 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 5 Mar 2019 13:03:02 +0000 Subject: [PATCH 23/33] Do convergence on 5 periods --- advection/weno.py | 117 ++++++---------------------------------------- 1 file changed, 13 insertions(+), 104 deletions(-) diff --git a/advection/weno.py b/advection/weno.py index 4f66faf..e5a1a81 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -287,95 +287,16 @@ def rk_substep_scipy(t, y): label="WENO3") -# #------------------------------------------------------------------------- -# # convergence test -# # Note that WENO schemes with standard weights lose convergence at -# # critical points. For high degree critical points they lose more orders. -# # The suggestion in Gerolymos is that you may expect to drop down to -# # order r-1 in the limit. -# # The Gaussian has all odd derivatives vanishing at the origin, so -# # the higher order schemes will lose accuracy. -# # For the Gaussian: -# # This shows clean 5th order convergence for r=3 -# # But for r=4-6 the best you get is ~6th order, and 5th order is more -# # realistic -# # For sin(x - sin(x)) type data Gerolymos expects better results -# # But the problem actually appears to be the time integrator -# # Switching to Dormand-Price 8th order from scipy (a hack) will make it -# # work for all cases. With sin(.. sin) data you get 2r - 2 thanks to -# # the one critical point. -# -# problem = "sine_sine" -# -# xmin =-1.0 -# xmax = 1.0 -## orders = [4] -# orders = [3, 4, 5, 6] -## N1 = [2**4*3**i//2**i for i in range(5)] -## N2 = [2**5*3**i//2**i for i in range(6)] -## N3 = [3**4*4**i//3**i for i in range(5)] -## N4 = [2**(4+i) for i in range(4)] -## N = numpy.unique(numpy.array(N1+N2+N3+N4, dtype=numpy.int)) -## N.sort() -## N = [32, 64, 128, 256, 512] -## N = [32, 64, 128] -# N = [24, 32, 54, 64, 81, 108, 128] -# -# errs = [] -# errsM = [] -# -# u = 1.0 -# -# colors="bygrc" -# -# for order in orders: -# ng = order+1 -# errs.append([]) -# errsM.append([]) -# for nx in N: -# print(order, nx) -# gu = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) -# su = WENOSimulation(gu, u, C=0.5, weno_order=order) -## guM = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) -## suM = WENOMSimulation(guM, u, C=0.5, weno_order=order) -# -# su.init_cond("sine_sine") -## suM.init_cond("sine_sine") -# ainit = su.grid.a.copy() -# -# su.evolve_scipy(num_periods=1) -## suM.evolve_scipy(num_periods=1) -# -# errs[-1].append(gu.norm(gu.a - ainit)) -## errsM[-1].append(guM.norm(guM.a - ainit)) -# -# pyplot.clf() -# N = numpy.array(N, dtype=numpy.float64) -# for n_order, order in enumerate(orders): -# pyplot.scatter(N, errs[n_order], -# color=colors[n_order], -# label=r"WENO, $r={}$".format(order)) -## pyplot.scatter(N, errsM[n_order], -## color=colors[n_order], -## label=r"WENOM, $r={}$".format(order)) -# pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-2), -# linestyle="--", color=colors[n_order], -# label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-2)) -## pyplot.plot(N, errs[n_order][len(N)-1]*(N[len(N)-1]/N)**4, -## color="k", label=r"$\mathcal{O}(\Delta x^4)$") -# -# ax = pyplot.gca() -# ax.set_ylim(numpy.min(errs)/5, numpy.max(errs)*5) -# ax.set_xscale('log') -# ax.set_yscale('log') -# -# pyplot.xlabel("N") -# pyplot.ylabel(r"$\| a^\mathrm{final} - a^\mathrm{init} \|_2$", -# fontsize=16) -# -# pyplot.legend(frameon=False) -# pyplot.savefig("weno-converge-sine-sine.pdf") -## pyplot.show() + #------------------------------------------------------------------------- + # convergence test + # Note that WENO schemes with standard weights lose convergence at + # critical points. For high degree critical points they lose more orders. + # The suggestion in Gerolymos is that you may expect to drop down to + # order r-1 in the limit. + # + # For the odd r values and using sine initial data we can get optimal + # convergence using 8th order time integration. For other cases the + # results are not so nice. #-------------- RK4 @@ -431,7 +352,7 @@ def rk_substep_scipy(t, y): pyplot.legend(frameon=False) pyplot.savefig("weno-converge-gaussian-rk4.pdf") -# pyplot.show() + pyplot.show() #-------------- Sine wave, 8th order time integrator @@ -441,7 +362,6 @@ def rk_substep_scipy(t, y): xmax = 1.0 orders = [3, 5, 7] N = [24, 32, 54, 64, 81, 108, 128] -# N = [32, 64, 108, 128] errs = [] errsM = [] @@ -458,18 +378,12 @@ def rk_substep_scipy(t, y): print(order, nx) gu = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) su = WENOSimulation(gu, u, C=0.5, weno_order=order) -# guM = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) -# suM = WENOMSimulation(guM, u, C=0.5, weno_order=order) su.init_cond("sine") -# suM.init_cond("sine") ainit = su.grid.a.copy() - su.evolve_scipy(num_periods=1) -# suM.evolve_scipy(num_periods=1) - + su.evolve_scipy(num_periods=5) errs[-1].append(gu.norm(gu.a - ainit, norm=2)) -# errsM[-1].append(guM.norm(guM.a - ainit, norm=2)) pyplot.clf() N = numpy.array(N, dtype=numpy.float64) @@ -477,14 +391,9 @@ def rk_substep_scipy(t, y): pyplot.scatter(N, errs[n_order], color=colors[n_order], label=r"WENO, $r={}$".format(order)) -# pyplot.scatter(N, errsM[n_order], -# color=colors[n_order], -# label=r"WENOM, $r={}$".format(order)) pyplot.plot(N, errs[n_order][0]*(N[0]/N)**(2*order-1), linestyle="--", color=colors[n_order], label=r"$\mathcal{{O}}(\Delta x^{{{}}})$".format(2*order-1)) -# pyplot.plot(N, errs[n_order][len(N)-1]*(N[len(N)-1]/N)**4, -# color="k", label=r"$\mathcal{O}(\Delta x^4)$") ax = pyplot.gca() ax.set_ylim(numpy.min(errs)/5, numpy.max(errs)*5) @@ -499,5 +408,5 @@ def rk_substep_scipy(t, y): lgd = ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) pyplot.savefig("weno-converge-sine.pdf", bbox_extra_artists=(lgd,), bbox_inches='tight') -# pyplot.show() + pyplot.show() \ No newline at end of file From 20c30bcaf11ae459afc6e24aa8e18cec89c5b45e Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 5 Mar 2019 13:37:17 +0000 Subject: [PATCH 24/33] Fix legend --- advection/compare_schemes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/advection/compare_schemes.py b/advection/compare_schemes.py index f1c379b..908298f 100644 --- a/advection/compare_schemes.py +++ b/advection/compare_schemes.py @@ -98,7 +98,7 @@ def errs_dg(m=3, nx=8, limiter=None): fig, ax = pyplot.subplots(1, 1) for i_o, order in enumerate(weno_orders): ax.loglog(weno_times[i_o, :], weno_errs[i_o, :], f'{colors[i_o]}o-', - label=f'WENO, order={order}') + label=f'WENO, r={order}') colors = "gyc" for i_m, m in enumerate(dg_ms): ax.loglog(dg_nolimit_times[i_m, :], dg_nolimit_errs[i_m, :], From ef13447572e4af71720f51aae5e9a7f504de64b6 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 5 Mar 2019 15:45:32 +0000 Subject: [PATCH 25/33] Using numba to speed up the WENO code --- advection/compare_schemes.py | 54 +++++++++----- advection/weno.py | 132 +++++++++++++++++++++++++++++++++++ 2 files changed, 168 insertions(+), 18 deletions(-) diff --git a/advection/compare_schemes.py b/advection/compare_schemes.py index 908298f..334bc3c 100644 --- a/advection/compare_schemes.py +++ b/advection/compare_schemes.py @@ -7,18 +7,21 @@ """ import timeit -from memory_profiler import memory_usage +#from memory_profiler import memory_usage import numpy from matplotlib import pyplot import advection import weno import dg -def run_weno(order=2, nx=32): +def run_weno(order=2, nx=32, optimized=False): g = advection.Grid1d(nx, ng=order+1, xmin=0, xmax=1) s = weno.WENOSimulation(g, 1, C=0.5, weno_order=order) s.init_cond("sine") - s.evolve_scipy() + if optimized: + s.evolve_scipy_jit() + else: + s.evolve_scipy() def run_dg(m=3, nx=8, limiter=None): g = dg.Grid1d(nx, ng=1, xmin=0, xmax=1, m=m) @@ -26,26 +29,29 @@ def run_dg(m=3, nx=8, limiter=None): s.init_cond("sine") s.evolve_scipy() -def weno_string(order=2, nx=32): - return f"run_weno(order={order}, nx={nx})" +def weno_string(order=2, nx=32, optimized=False): + return f"run_weno(order={order}, nx={nx}, optimized={optimized})" def dg_string(m=3, nx=8, limiter=None): return f"run_dg(m={m}, nx={nx}, limiter='{limiter}')" -def time_weno(order=2, nx=32): - return timeit.timeit(weno_string(order, nx), +def time_weno(order=2, nx=32, optimized=False): + return timeit.timeit(weno_string(order, nx, optimized), globals=globals(), number=1) def time_dg(m=3, nx=8, limiter=None): return timeit.timeit(dg_string(m, nx, limiter), globals=globals(), number=1) -def errs_weno(order=2, nx=32): +def errs_weno(order=2, nx=32, optimized=False): g = advection.Grid1d(nx, ng=order+1, xmin=0, xmax=1) s = weno.WENOSimulation(g, 1, C=0.5, weno_order=order) s.init_cond("sine") a_init = g.a.copy() - s.evolve_scipy() + if optimized: + s.evolve_scipy_jit() + else: + s.evolve_scipy() return g.norm(g.a - a_init) def errs_dg(m=3, nx=8, limiter=None): @@ -61,15 +67,25 @@ def errs_dg(m=3, nx=8, limiter=None): #weno_orders = [3] #weno_N = [24, 32, 54] weno_times = numpy.zeros((len(weno_orders), len(weno_N))) +weno_opt_times = numpy.zeros((len(weno_orders), len(weno_N))) weno_errs = numpy.zeros_like(weno_times) -weno_memory = numpy.zeros_like(weno_times) +weno_opt_errs = numpy.zeros_like(weno_opt_times) +#weno_memory = numpy.zeros_like(weno_times) +#weno_opt_memory = numpy.zeros_like(weno_opt_times) + +# Do one evolution to kick-start numba +run_weno(3, 8, optimized=True) for i_o, order in enumerate(weno_orders): for i_n, nx in enumerate(weno_N): print(f"WENO{order}, {nx} points") weno_errs[i_o, i_n] = errs_weno(order, nx) weno_times[i_o, i_n] = time_weno(order, nx) - weno_memory[i_o, i_n] = max(memory_usage((run_weno, (order, nx)))) +# weno_memory[i_o, i_n] = max(memory_usage((run_weno, (order, nx)))) + weno_opt_errs[i_o, i_n] = errs_weno(order, nx, optimized=True) + weno_opt_times[i_o, i_n] = time_weno(order, nx, optimized=True) +# weno_opt_memory[i_o, i_n] = max(memory_usage((run_weno, +# (order, nx, optimized)))) dg_ms = [2, 3, 4] dg_N = 2**numpy.array(range(3, 9)) @@ -79,30 +95,32 @@ def errs_dg(m=3, nx=8, limiter=None): dg_moment_times = numpy.zeros((len(dg_ms), len(dg_N))) dg_nolimit_errs = numpy.zeros_like(dg_nolimit_times) dg_moment_errs = numpy.zeros_like(dg_moment_times) -dg_nolimit_memory = numpy.zeros_like(dg_nolimit_times) -dg_moment_memory = numpy.zeros_like(dg_moment_times) +#dg_nolimit_memory = numpy.zeros_like(dg_nolimit_times) +#dg_moment_memory = numpy.zeros_like(dg_moment_times) for i_m, m in enumerate(dg_ms): for i_n, nx in enumerate(dg_N): print(f"DG, m={m}, {nx} points") dg_nolimit_errs[i_m, i_n] = errs_dg(m, nx, None) dg_nolimit_times[i_m, i_n] = time_dg(m, nx, None) - dg_nolimit_memory[i_m, i_n] = max(memory_usage((run_dg, - (m, nx, None)))) +# dg_nolimit_memory[i_m, i_n] = max(memory_usage((run_dg, +# (m, nx, None)))) dg_moment_errs[i_m, i_n] = errs_dg(m, nx, 'moment') dg_moment_times[i_m, i_n] = time_dg(m, nx, 'moment') - dg_moment_memory[i_m, i_n] = max(memory_usage((run_dg, - (m, nx, 'moment')))) +# dg_moment_memory[i_m, i_n] = max(memory_usage((run_dg, +# (m, nx, 'moment')))) colors = "brk" fig, ax = pyplot.subplots(1, 1) for i_o, order in enumerate(weno_orders): ax.loglog(weno_times[i_o, :], weno_errs[i_o, :], f'{colors[i_o]}o-', label=f'WENO, r={order}') + ax.loglog(weno_opt_times[i_o, :], weno_opt_errs[i_o, :], + f'{colors[i_o]}o--', label=f'WENO, optimized, r={order}') colors = "gyc" for i_m, m in enumerate(dg_ms): ax.loglog(dg_nolimit_times[i_m, :], dg_nolimit_errs[i_m, :], - f'{colors[i_m]}s--', label=rf'DG, no limiter, $m={m}$') + f'{colors[i_m]}s-.', label=rf'DG, no limiter, $m={m}$') ax.loglog(dg_moment_times[i_m, :], dg_moment_errs[i_m, :], f'{colors[i_m]}^:', label=rf'DG, moment limiter, $m={m}$') ax.set_xlabel("Runtime [s]") diff --git a/advection/weno.py b/advection/weno.py index c73550b..125cae1 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -4,6 +4,8 @@ import weno_coefficients from scipy.integrate import ode +from numba import jit + def weno(order, q): """ @@ -47,6 +49,49 @@ def weno(order, q): return qL +@jit +def weno_jit(order, q): + """ + Do WENO reconstruction + + Parameters + ---------- + + order : int + The stencil width + q : numpy array + Scalar data to reconstruct + + Returns + ------- + + qL : numpy array + Reconstructed data - boundary points are zero + """ + C = weno_coefficients.C_all[order] + a = weno_coefficients.a_all[order] + sigma = weno_coefficients.sigma_all[order] + + qL = numpy.zeros_like(q) + beta = numpy.zeros((order, len(q))) + w = numpy.zeros_like(beta) + np = len(q) - 2 * order + epsilon = 1e-16 + for i in range(order, np+order): + q_stencils = numpy.zeros(order) + alpha = numpy.zeros(order) + for k in range(order): + for l in range(order): + for m in range(l+1): + beta[k, i] += sigma[k, l, m] * q[i+k-l] * q[i+k-m] + alpha[k] = C[k] / (epsilon + beta[k, i]**2) + for l in range(order): + q_stencils[k] += a[k, l] * q[i+k-l] + w[:, i] = alpha / numpy.sum(alpha) + qL[i] = numpy.dot(w[:, i], q_stencils) + + return qL + def weno_M(order, q): """ @@ -131,6 +176,25 @@ def rk_substep(self): rhs[1:-1] = 1/g.dx * (flux[1:-1] - flux[2:]) return rhs + @jit + def rk_substep_jit(self): + + g = self.grid + g.fill_BCs() + f = self.u * g.a + alpha = abs(self.u) + fp = (f + alpha * g.a) / 2 + fm = (f - alpha * g.a) / 2 + fpr = g.scratch_array() + fml = g.scratch_array() + flux = g.scratch_array() + fpr[1:] = weno_jit(self.weno_order, fp[:-1]) + fml[-1::-1] = weno_jit(self.weno_order, fm[-1::-1]) + flux[1:-1] = fpr[1:-1] + fml[1:-1] + rhs = g.scratch_array() + rhs[1:-1] = 1/g.dx * (flux[1:-1] - flux[2:]) + return rhs + def evolve(self, num_periods=1): """ evolve the linear advection equation using RK4 """ @@ -165,6 +229,40 @@ def evolve(self, num_periods=1): self.t += dt + @jit + def evolve_jit(self, num_periods=1): + """ evolve the linear advection equation using RK4 """ + self.t = 0.0 + g = self.grid + + tmax = num_periods*self.period() + + # main evolution loop + while self.t < tmax: + + # fill the boundary conditions + g.fill_BCs() + + # get the timestep + dt = self.timestep() + + if self.t + dt > tmax: + dt = tmax - self.t + + # RK4 + # Store the data at the start of the step + a_start = g.a.copy() + k1 = dt * self.rk_substep_jit() + g.a = a_start + k1 / 2 + k2 = dt * self.rk_substep_jit() + g.a = a_start + k2 / 2 + k3 = dt * self.rk_substep_jit() + g.a = a_start + k3 + k4 = dt * self.rk_substep_jit() + g.a = a_start + (k1 + 2 * (k2 + k3) + k4) / 6 + + self.t += dt + def evolve_scipy(self, num_periods=1): """ evolve the linear advection equation using RK4 """ @@ -200,6 +298,40 @@ def rk_substep_scipy(t, y): r.integrate(r.t+dt) g.a[:] = r.y + def evolve_scipy_jit(self, num_periods=1): + """ evolve the linear advection equation using RK4 """ + self.t = 0.0 + g = self.grid + + def rk_substep_scipy_jit(t, y): + # Periodic BCs + y[:g.ng] = y[-2*g.ng:-g.ng] + y[-g.ng:] = y[g.ng:2*g.ng] + f = self.u * y + alpha = abs(self.u) + fp = (f + alpha * y) / 2 + fm = (f - alpha * y) / 2 + fpr = g.scratch_array() + fml = g.scratch_array() + flux = g.scratch_array() + fpr[1:] = weno_jit(self.weno_order, fp[:-1]) + fml[-1::-1] = weno_jit(self.weno_order, fm[-1::-1]) + flux[1:-1] = fpr[1:-1] + fml[1:-1] + rhs = g.scratch_array() + rhs[1:-1] = 1/g.dx * (flux[1:-1] - flux[2:]) + return rhs + + tmax = num_periods*self.period() + r = ode(rk_substep_scipy_jit).set_integrator('dop853') + r.set_initial_value(g.a, 0) + dt = self.timestep() + + # main evolution loop + while r.successful() and r.t < tmax: + dt = min(dt, tmax - r.t) + r.integrate(r.t+dt) + g.a[:] = r.y + class WENOMSimulation(WENOSimulation): From dbaba5aae7be9abf549aa484f7b2c56e86107a88 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 6 Mar 2019 09:45:26 +0000 Subject: [PATCH 26/33] Redo efficiency with only limited DG and everything optimized --- advection/compare_schemes.py | 84 +++++++--------- advection/dg.py | 187 +++++++++++++++++++++++++++++++++++ advection/weno.py | 185 +++++----------------------------- 3 files changed, 249 insertions(+), 207 deletions(-) diff --git a/advection/compare_schemes.py b/advection/compare_schemes.py index 334bc3c..2250bea 100644 --- a/advection/compare_schemes.py +++ b/advection/compare_schemes.py @@ -7,74 +7,76 @@ """ import timeit -#from memory_profiler import memory_usage +# from memory_profiler import memory_usage import numpy from matplotlib import pyplot import advection import weno import dg -def run_weno(order=2, nx=32, optimized=False): + +def run_weno(order=2, nx=32): g = advection.Grid1d(nx, ng=order+1, xmin=0, xmax=1) s = weno.WENOSimulation(g, 1, C=0.5, weno_order=order) s.init_cond("sine") - if optimized: - s.evolve_scipy_jit() - else: - s.evolve_scipy() - + s.evolve_scipy() + + def run_dg(m=3, nx=8, limiter=None): g = dg.Grid1d(nx, ng=1, xmin=0, xmax=1, m=m) s = dg.Simulation(g, 1, C=0.5/(2*m+1), limiter=limiter) s.init_cond("sine") s.evolve_scipy() -def weno_string(order=2, nx=32, optimized=False): - return f"run_weno(order={order}, nx={nx}, optimized={optimized})" + +def weno_string(order=2, nx=32): + return f"run_weno(order={order}, nx={nx})" + def dg_string(m=3, nx=8, limiter=None): return f"run_dg(m={m}, nx={nx}, limiter='{limiter}')" -def time_weno(order=2, nx=32, optimized=False): - return timeit.timeit(weno_string(order, nx, optimized), + +def time_weno(order=2, nx=32): + return timeit.timeit(weno_string(order, nx), globals=globals(), number=1) + def time_dg(m=3, nx=8, limiter=None): return timeit.timeit(dg_string(m, nx, limiter), globals=globals(), number=1) -def errs_weno(order=2, nx=32, optimized=False): + +def errs_weno(order=2, nx=32): g = advection.Grid1d(nx, ng=order+1, xmin=0, xmax=1) s = weno.WENOSimulation(g, 1, C=0.5, weno_order=order) s.init_cond("sine") a_init = g.a.copy() - if optimized: - s.evolve_scipy_jit() - else: - s.evolve_scipy() + s.evolve_scipy() return g.norm(g.a - a_init) - + + def errs_dg(m=3, nx=8, limiter=None): g = dg.Grid1d(nx, ng=1, xmin=0, xmax=1, m=m) s = dg.Simulation(g, 1, C=0.5/(2*m+1), limiter=limiter) s.init_cond("sine") a_init = g.a.copy() - s.evolve_scipy() + s.evolve_scipy_jit() return g.norm(g.a - a_init) + weno_orders = [3, 5, 7] -weno_N = [8, 12, 16, 24, 32, 54, 64] -#weno_orders = [3] -#weno_N = [24, 32, 54] +weno_N = [12, 16, 24, 32, 54, 64, 96, 128] +# weno_orders = [3] +# weno_N = [24, 32, 54] weno_times = numpy.zeros((len(weno_orders), len(weno_N))) -weno_opt_times = numpy.zeros((len(weno_orders), len(weno_N))) weno_errs = numpy.zeros_like(weno_times) -weno_opt_errs = numpy.zeros_like(weno_opt_times) -#weno_memory = numpy.zeros_like(weno_times) -#weno_opt_memory = numpy.zeros_like(weno_opt_times) +# weno_memory = numpy.zeros_like(weno_times) +# weno_opt_memory = numpy.zeros_like(weno_opt_times) # Do one evolution to kick-start numba -run_weno(3, 8, optimized=True) +run_weno(3, 8) +run_dg(3, 8) for i_o, order in enumerate(weno_orders): for i_n, nx in enumerate(weno_N): @@ -82,29 +84,19 @@ def errs_dg(m=3, nx=8, limiter=None): weno_errs[i_o, i_n] = errs_weno(order, nx) weno_times[i_o, i_n] = time_weno(order, nx) # weno_memory[i_o, i_n] = max(memory_usage((run_weno, (order, nx)))) - weno_opt_errs[i_o, i_n] = errs_weno(order, nx, optimized=True) - weno_opt_times[i_o, i_n] = time_weno(order, nx, optimized=True) -# weno_opt_memory[i_o, i_n] = max(memory_usage((run_weno, -# (order, nx, optimized)))) - -dg_ms = [2, 3, 4] -dg_N = 2**numpy.array(range(3, 9)) -#dg_ms = [2] -#dg_N = 2**numpy.array(range(3, 6)) -dg_nolimit_times = numpy.zeros((len(dg_ms), len(dg_N))) + +dg_ms = [2, 4, 8] +dg_N = 2**numpy.array(range(3, 7)) +# dg_ms = [2] +# dg_N = 2**numpy.array(range(3, 6)) dg_moment_times = numpy.zeros((len(dg_ms), len(dg_N))) -dg_nolimit_errs = numpy.zeros_like(dg_nolimit_times) dg_moment_errs = numpy.zeros_like(dg_moment_times) -#dg_nolimit_memory = numpy.zeros_like(dg_nolimit_times) -#dg_moment_memory = numpy.zeros_like(dg_moment_times) +# dg_nolimit_memory = numpy.zeros_like(dg_nolimit_times) +# dg_moment_memory = numpy.zeros_like(dg_moment_times) for i_m, m in enumerate(dg_ms): for i_n, nx in enumerate(dg_N): print(f"DG, m={m}, {nx} points") - dg_nolimit_errs[i_m, i_n] = errs_dg(m, nx, None) - dg_nolimit_times[i_m, i_n] = time_dg(m, nx, None) -# dg_nolimit_memory[i_m, i_n] = max(memory_usage((run_dg, -# (m, nx, None)))) dg_moment_errs[i_m, i_n] = errs_dg(m, nx, 'moment') dg_moment_times[i_m, i_n] = time_dg(m, nx, 'moment') # dg_moment_memory[i_m, i_n] = max(memory_usage((run_dg, @@ -115,14 +107,10 @@ def errs_dg(m=3, nx=8, limiter=None): for i_o, order in enumerate(weno_orders): ax.loglog(weno_times[i_o, :], weno_errs[i_o, :], f'{colors[i_o]}o-', label=f'WENO, r={order}') - ax.loglog(weno_opt_times[i_o, :], weno_opt_errs[i_o, :], - f'{colors[i_o]}o--', label=f'WENO, optimized, r={order}') colors = "gyc" for i_m, m in enumerate(dg_ms): - ax.loglog(dg_nolimit_times[i_m, :], dg_nolimit_errs[i_m, :], - f'{colors[i_m]}s-.', label=rf'DG, no limiter, $m={m}$') ax.loglog(dg_moment_times[i_m, :], dg_moment_errs[i_m, :], - f'{colors[i_m]}^:', label=rf'DG, moment limiter, $m={m}$') + f'{colors[i_m]}^:', label=rf'DG, $m={m}$') ax.set_xlabel("Runtime [s]") ax.set_ylabel(r"$\|$Error$\|_2$") ax.set_title("Efficiency of WENO vs DG") diff --git a/advection/dg.py b/advection/dg.py index 6b02ca7..327de7b 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -8,6 +8,7 @@ import matplotlib as mpl import quadpy from scipy.integrate import ode +from numba import jit mpl.rcParams['mathtext.fontset'] = 'cm' mpl.rcParams['mathtext.rm'] = 'serif' @@ -148,6 +149,141 @@ def minmod3(a1, a2, a3): return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) +@jit +def minmod3_jit(a1, a2, a3): + """ + Utility function that does minmod on three inputs + """ + signs1 = a1 * a2 + signs2 = a1 * a3 + signs3 = a2 * a3 + same_sign = numpy.logical_and(numpy.logical_and(signs1 > 0, + signs2 > 0), + signs3 > 0) + minmod = numpy.min(numpy.abs(numpy.vstack((a1, a2, a3))), + axis=0) * numpy.sign(a1) + + return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) + + +def limit(g, limiter): + """ + After evolution, limit the slopes. + """ + + # Limiting! + + if limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[1:, :] = 0 + a_cell = (g.V @ a_zeromode)[0, :] + a_minus = g.a[0, :] + a_plus = g.a[-1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[i+1, 1:-1] + a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) + a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor + did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], + updated_mode[1:-1]) + a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + numpy.logical_not(did_it_limit)) + # Get back nodal values + limited_a = g.modal_to_nodal(a_modal) + + else: + limited_a = g.a + + return limited_a + + +@jit +def limit_jit(g, limiter): + """ + After evolution, limit the slopes. + """ + + # Limiting! + + if limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[1:, :] = 0 + a_cell = (g.V @ a_zeromode)[0, :] + a_minus = g.a[0, :] + a_plus = g.a[-1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[i+1, 1:-1] + a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) + a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor + did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], + updated_mode[1:-1]) + a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + numpy.logical_not(did_it_limit)) + # Get back nodal values + limited_a = g.modal_to_nodal(a_modal) + + return limited_a + + class Simulation(object): def __init__(self, grid, u, C=0.8, limiter=None): @@ -386,6 +522,57 @@ def rk_substep_scipy(t, y): g.a = self.limit() self.t += dt + def evolve_scipy_jit(self, num_periods=1): + """ evolve the linear advection equation using scipy """ + self.t = 0.0 + g = self.grid + + def rk_substep_scipy_jit(t, y): + local_a = numpy.reshape(y, g.a.shape) + # Periodic BCs + local_a[:, :g.ng] = local_a[:, -2*g.ng:-g.ng] + local_a[:, -g.ng:] = local_a[:, g.ng:2*g.ng] + # Integrate flux over element + f = self.u * local_a + interior_f = g.S.T @ f + # Use Riemann solver to get fluxes between elements + al = numpy.zeros(g.nx+2*g.ng) + ar = numpy.zeros(g.nx+2*g.ng) + # i is looping over interfaces, so al is the right edge of the left + # element, etc. + for i in range(g.ilo, g.ihi+2): + al[i] = local_a[-1, i-1] + ar[i] = local_a[0, i] + boundary_f = self.riemann(al, ar) + rhs = interior_f + rhs[0, 1:-1] += boundary_f[1:-1] + rhs[-1, 1:-1] -= boundary_f[2:] + + # Multiply by mass matrix (inverse). + rhs_i = 2 / g.dx * g.M_inv @ rhs + + return numpy.ravel(rhs_i, order='C') + + tmax = num_periods*self.period() + + # main evolution loop + while self.t < tmax: + # fill the boundary conditions + g.fill_BCs() + + # get the timestep + dt = self.timestep() + + if self.t + dt > tmax: + dt = tmax - self.t + + r = ode(rk_substep_scipy_jit).set_integrator('dop853') + r.set_initial_value(numpy.ravel(g.a), self.t) + r.integrate(r.t+dt) + g.a[:, :] = numpy.reshape(r.y, g.a.shape) + g.a = limit(g, self.limiter) + self.t += dt + if __name__ == "__main__": diff --git a/advection/weno.py b/advection/weno.py index 125cae1..504bc59 100644 --- a/advection/weno.py +++ b/advection/weno.py @@ -7,64 +7,22 @@ from numba import jit +@jit def weno(order, q): """ Do WENO reconstruction - - Parameters - ---------- - - order : int - The stencil width - q : numpy array - Scalar data to reconstruct - - Returns - ------- - - qL : numpy array - Reconstructed data - boundary points are zero - """ - C = weno_coefficients.C_all[order] - a = weno_coefficients.a_all[order] - sigma = weno_coefficients.sigma_all[order] - - qL = numpy.zeros_like(q) - beta = numpy.zeros((order, len(q))) - w = numpy.zeros_like(beta) - np = len(q) - 2 * order - epsilon = 1e-16 - for i in range(order, np+order): - q_stencils = numpy.zeros(order) - alpha = numpy.zeros(order) - for k in range(order): - for l in range(order): - for m in range(l+1): - beta[k, i] += sigma[k, l, m] * q[i+k-l] * q[i+k-m] - alpha[k] = C[k] / (epsilon + beta[k, i]**2) - for l in range(order): - q_stencils[k] += a[k, l] * q[i+k-l] - w[:, i] = alpha / numpy.sum(alpha) - qL[i] = numpy.dot(w[:, i], q_stencils) - - return qL -@jit -def weno_jit(order, q): - """ - Do WENO reconstruction - Parameters ---------- - + order : int The stencil width q : numpy array Scalar data to reconstruct - + Returns ------- - + qL : numpy array Reconstructed data - boundary points are zero """ @@ -89,25 +47,25 @@ def weno_jit(order, q): q_stencils[k] += a[k, l] * q[i+k-l] w[:, i] = alpha / numpy.sum(alpha) qL[i] = numpy.dot(w[:, i], q_stencils) - + return qL def weno_M(order, q): """ Do WENOM reconstruction following Gerolymos equation (18) - + Parameters ---------- - + order : int The stencil width q : numpy array Scalar data to reconstruct - + Returns ------- - + qL : numpy array Reconstructed data - boundary points are zero """ @@ -135,29 +93,28 @@ def weno_M(order, q): (C**2 + w_JS * (1 - 2 * C)) w[:, i] = alpha / numpy.sum(alpha) qL[i] = numpy.dot(w[:, i], q_stencils) - + return qL class WENOSimulation(advection.Simulation): - + def __init__(self, grid, u, C=0.8, weno_order=3): self.grid = grid - self.t = 0.0 # simulation time - self.u = u # the constant advective velocity - self.C = C # CFL number + self.t = 0.0 # simulation time + self.u = u # the constant advective velocity + self.C = C # CFL number self.weno_order = weno_order - def init_cond(self, type="tophat"): """ initialize the data """ if type == "sine_sine": - self.grid.a[:] = numpy.sin(numpy.pi*self.grid.x - + self.grid.a[:] = numpy.sin(numpy.pi*self.grid.x - numpy.sin(numpy.pi*self.grid.x) / numpy.pi) else: super().init_cond(type) - + @jit def rk_substep(self): g = self.grid @@ -177,25 +134,6 @@ def rk_substep(self): return rhs @jit - def rk_substep_jit(self): - - g = self.grid - g.fill_BCs() - f = self.u * g.a - alpha = abs(self.u) - fp = (f + alpha * g.a) / 2 - fm = (f - alpha * g.a) / 2 - fpr = g.scratch_array() - fml = g.scratch_array() - flux = g.scratch_array() - fpr[1:] = weno_jit(self.weno_order, fp[:-1]) - fml[-1::-1] = weno_jit(self.weno_order, fm[-1::-1]) - flux[1:-1] = fpr[1:-1] + fml[1:-1] - rhs = g.scratch_array() - rhs[1:-1] = 1/g.dx * (flux[1:-1] - flux[2:]) - return rhs - - def evolve(self, num_periods=1): """ evolve the linear advection equation using RK4 """ self.t = 0.0 @@ -229,46 +167,11 @@ def evolve(self, num_periods=1): self.t += dt - @jit - def evolve_jit(self, num_periods=1): - """ evolve the linear advection equation using RK4 """ - self.t = 0.0 - g = self.grid - - tmax = num_periods*self.period() - - # main evolution loop - while self.t < tmax: - - # fill the boundary conditions - g.fill_BCs() - - # get the timestep - dt = self.timestep() - - if self.t + dt > tmax: - dt = tmax - self.t - - # RK4 - # Store the data at the start of the step - a_start = g.a.copy() - k1 = dt * self.rk_substep_jit() - g.a = a_start + k1 / 2 - k2 = dt * self.rk_substep_jit() - g.a = a_start + k2 / 2 - k3 = dt * self.rk_substep_jit() - g.a = a_start + k3 - k4 = dt * self.rk_substep_jit() - g.a = a_start + (k1 + 2 * (k2 + k3) + k4) / 6 - - self.t += dt - - def evolve_scipy(self, num_periods=1): """ evolve the linear advection equation using RK4 """ self.t = 0.0 g = self.grid - + def rk_substep_scipy(t, y): # Periodic BCs y[:g.ng] = y[-2*g.ng:-g.ng] @@ -298,45 +201,11 @@ def rk_substep_scipy(t, y): r.integrate(r.t+dt) g.a[:] = r.y - def evolve_scipy_jit(self, num_periods=1): - """ evolve the linear advection equation using RK4 """ - self.t = 0.0 - g = self.grid - - def rk_substep_scipy_jit(t, y): - # Periodic BCs - y[:g.ng] = y[-2*g.ng:-g.ng] - y[-g.ng:] = y[g.ng:2*g.ng] - f = self.u * y - alpha = abs(self.u) - fp = (f + alpha * y) / 2 - fm = (f - alpha * y) / 2 - fpr = g.scratch_array() - fml = g.scratch_array() - flux = g.scratch_array() - fpr[1:] = weno_jit(self.weno_order, fp[:-1]) - fml[-1::-1] = weno_jit(self.weno_order, fm[-1::-1]) - flux[1:-1] = fpr[1:-1] + fml[1:-1] - rhs = g.scratch_array() - rhs[1:-1] = 1/g.dx * (flux[1:-1] - flux[2:]) - return rhs - - tmax = num_periods*self.period() - r = ode(rk_substep_scipy_jit).set_integrator('dop853') - r.set_initial_value(g.a, 0) - dt = self.timestep() - - # main evolution loop - while r.successful() and r.t < tmax: - dt = min(dt, tmax - r.t) - r.integrate(r.t+dt) - g.a[:] = r.y - class WENOMSimulation(WENOSimulation): def rk_substep(self): - + g = self.grid g.fill_BCs() f = self.u * g.a @@ -352,13 +221,12 @@ def rk_substep(self): rhs = g.scratch_array() rhs[1:-1] = 1/g.dx * (flux[1:-1] - flux[2:]) return rhs - - + def evolve_scipy(self, num_periods=1): """ evolve the linear advection equation using scipy """ self.t = 0.0 g = self.grid - + def rk_substep_scipy(t, y): # Periodic BCs y[:g.ng] = y[-2*g.ng:-g.ng] @@ -387,12 +255,12 @@ def rk_substep_scipy(t, y): dt = min(dt, tmax - r.t) r.integrate(r.t+dt) g.a[:] = r.y - + if __name__ == "__main__": - #------------------------------------------------------------------------- + # ------------------------------------------------------------------------- # compute WENO3 case xmin = 0.0 @@ -404,7 +272,7 @@ def rk_substep_scipy(t, y): g = advection.Grid1d(nx, ng, xmin=xmin, xmax=xmax) u = 1.0 - + s = WENOSimulation(g, u, C=0.5, weno_order=3) s.init_cond("gaussian") @@ -413,12 +281,11 @@ def rk_substep_scipy(t, y): s.evolve(num_periods=1) pyplot.plot(g.x[g.ilo:g.ihi+1], ainit[g.ilo:g.ihi+1], - ls=":", label="exact") + ls=":", label="exact") pyplot.plot(g.x[g.ilo:g.ihi+1], g.a[g.ilo:g.ihi+1], - label="WENO3") - - + label="WENO3") + # #------------------------------------------------------------------------- # # convergence test # # Note that WENO schemes with standard weights lose convergence at From c58bda8caab561b175fb54ceec482d0f9f0e8dd7 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 6 Mar 2019 10:42:50 +0000 Subject: [PATCH 27/33] Annotate plot and tidy DG optimization --- advection/compare_schemes.py | 18 ++++- advection/dg.py | 126 +---------------------------------- 2 files changed, 17 insertions(+), 127 deletions(-) diff --git a/advection/compare_schemes.py b/advection/compare_schemes.py index 2250bea..a14140d 100644 --- a/advection/compare_schemes.py +++ b/advection/compare_schemes.py @@ -61,7 +61,7 @@ def errs_dg(m=3, nx=8, limiter=None): s = dg.Simulation(g, 1, C=0.5/(2*m+1), limiter=limiter) s.init_cond("sine") a_init = g.a.copy() - s.evolve_scipy_jit() + s.evolve_scipy() return g.norm(g.a - a_init) @@ -72,7 +72,6 @@ def errs_dg(m=3, nx=8, limiter=None): weno_times = numpy.zeros((len(weno_orders), len(weno_N))) weno_errs = numpy.zeros_like(weno_times) # weno_memory = numpy.zeros_like(weno_times) -# weno_opt_memory = numpy.zeros_like(weno_opt_times) # Do one evolution to kick-start numba run_weno(3, 8) @@ -114,6 +113,21 @@ def errs_dg(m=3, nx=8, limiter=None): ax.set_xlabel("Runtime [s]") ax.set_ylabel(r"$\|$Error$\|_2$") ax.set_title("Efficiency of WENO vs DG") + +i_o = 0 +i_n = 4 +con_style = "arc,angleA=180,armA=50,rad=10" +ax.annotate(fr'WENO, $N={weno_N[i_n]}$', textcoords='data', + xy=(weno_times[i_o, i_n], 1.5*weno_errs[i_o, i_n]), + xytext=(5e-1, 1e-2), + arrowprops=dict(arrowstyle='->', connectionstyle=con_style)) +i_m = 1 +i_n = 0 +ax.annotate(fr'DG, $N={dg_N[i_n]}$', textcoords='data', + xy=(dg_moment_times[i_m, i_n], 0.5*dg_moment_errs[i_m, i_n]), + xytext=(3e-1, 1e-8), + arrowprops=dict(arrowstyle='->', connectionstyle=con_style)) + fig.tight_layout() lgd = ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) fig.savefig('dg_weno_efficiency.pdf', diff --git a/advection/dg.py b/advection/dg.py index 327de7b..16e0343 100644 --- a/advection/dg.py +++ b/advection/dg.py @@ -133,6 +133,7 @@ def norm(self, e): return numpy.sqrt(self.dx*numpy.sum(e[:, self.ilo:self.ihi+1]**2)) +@jit def minmod3(a1, a2, a3): """ Utility function that does minmod on three inputs @@ -150,22 +151,6 @@ def minmod3(a1, a2, a3): @jit -def minmod3_jit(a1, a2, a3): - """ - Utility function that does minmod on three inputs - """ - signs1 = a1 * a2 - signs2 = a1 * a3 - signs3 = a2 * a3 - same_sign = numpy.logical_and(numpy.logical_and(signs1 > 0, - signs2 > 0), - signs3 > 0) - minmod = numpy.min(numpy.abs(numpy.vstack((a1, a2, a3))), - axis=0) * numpy.sign(a1) - - return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) - - def limit(g, limiter): """ After evolution, limit the slopes. @@ -226,64 +211,6 @@ def limit(g, limiter): return limited_a -@jit -def limit_jit(g, limiter): - """ - After evolution, limit the slopes. - """ - - # Limiting! - - if limiter == "moment": - - # Limiting, using moment limiting (Hesthaven p 445-7) - theta = 2 - a_modal = g.nodal_to_modal() - # First, work out where limiting is needed - limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) - limiting_todo[:g.ng] = False - limiting_todo[-g.ng:] = False - # Get the cell average and the nodal values at the boundaries - a_zeromode = a_modal.copy() - a_zeromode[1:, :] = 0 - a_cell = (g.V @ a_zeromode)[0, :] - a_minus = g.a[0, :] - a_plus = g.a[-1, :] - # From the cell averages and boundary values we can construct - # alternate values at the boundaries - a_left = numpy.zeros(g.nx+2*g.ng) - a_right = numpy.zeros(g.nx+2*g.ng) - a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], - a_cell[1:-1] - a_cell[:-2], - a_cell[2:] - a_cell[1:-1]) - a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], - a_cell[1:-1] - a_cell[:-2], - a_cell[2:] - a_cell[1:-1]) - limiting_todo[1:-1] = numpy.logical_not( - numpy.logical_and(numpy.isclose(a_left[1:-1], - a_minus[1:-1]), - numpy.isclose(a_right[1:-1], - a_plus[1:-1]))) - # Now, do the limiting. Modify moments where needed, and as soon as - # limiting isn't needed, stop - updated_mode = numpy.zeros(g.nx+2*g.ng) - for i in range(g.m-1, -1, -1): - factor = numpy.sqrt((2*i+3)*(2*i+5)) - a1 = factor * a_modal[i+1, 1:-1] - a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) - a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) - updated_mode[1:-1] = minmod3(a1, a2, a3) / factor - did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], - updated_mode[1:-1]) - a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] - limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], - numpy.logical_not(did_it_limit)) - # Get back nodal values - limited_a = g.modal_to_nodal(a_modal) - - return limited_a - - class Simulation(object): def __init__(self, grid, u, C=0.8, limiter=None): @@ -522,57 +449,6 @@ def rk_substep_scipy(t, y): g.a = self.limit() self.t += dt - def evolve_scipy_jit(self, num_periods=1): - """ evolve the linear advection equation using scipy """ - self.t = 0.0 - g = self.grid - - def rk_substep_scipy_jit(t, y): - local_a = numpy.reshape(y, g.a.shape) - # Periodic BCs - local_a[:, :g.ng] = local_a[:, -2*g.ng:-g.ng] - local_a[:, -g.ng:] = local_a[:, g.ng:2*g.ng] - # Integrate flux over element - f = self.u * local_a - interior_f = g.S.T @ f - # Use Riemann solver to get fluxes between elements - al = numpy.zeros(g.nx+2*g.ng) - ar = numpy.zeros(g.nx+2*g.ng) - # i is looping over interfaces, so al is the right edge of the left - # element, etc. - for i in range(g.ilo, g.ihi+2): - al[i] = local_a[-1, i-1] - ar[i] = local_a[0, i] - boundary_f = self.riemann(al, ar) - rhs = interior_f - rhs[0, 1:-1] += boundary_f[1:-1] - rhs[-1, 1:-1] -= boundary_f[2:] - - # Multiply by mass matrix (inverse). - rhs_i = 2 / g.dx * g.M_inv @ rhs - - return numpy.ravel(rhs_i, order='C') - - tmax = num_periods*self.period() - - # main evolution loop - while self.t < tmax: - # fill the boundary conditions - g.fill_BCs() - - # get the timestep - dt = self.timestep() - - if self.t + dt > tmax: - dt = tmax - self.t - - r = ode(rk_substep_scipy_jit).set_integrator('dop853') - r.set_initial_value(numpy.ravel(g.a), self.t) - r.integrate(r.t+dt) - g.a[:, :] = numpy.reshape(r.y, g.a.shape) - g.a = limit(g, self.limiter) - self.t += dt - if __name__ == "__main__": From 24bd2bed5693baefc8f707d878e0a095e99640b4 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 6 Mar 2019 13:10:02 +0000 Subject: [PATCH 28/33] DG for Burgers. Cross check BCs in advection --- burgers/dg_burgers.py | 432 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 432 insertions(+) create mode 100644 burgers/dg_burgers.py diff --git a/burgers/dg_burgers.py b/burgers/dg_burgers.py new file mode 100644 index 0000000..e9a7bef --- /dev/null +++ b/burgers/dg_burgers.py @@ -0,0 +1,432 @@ +""" +Discontinuous Galerkin for Burgers equation. +""" + +import numpy +from numpy.polynomial import legendre +from matplotlib import pyplot +import matplotlib as mpl +import quadpy +from numba import jit +import tqdm + +mpl.rcParams['mathtext.fontset'] = 'cm' +mpl.rcParams['mathtext.rm'] = 'serif' + + +class Grid1d(object): + + def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): + + assert m > 0 + + self.ng = ng + self.nx = nx + self.m = m + + self.xmin = xmin + self.xmax = xmax + + # 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 + (numpy.arange(nx+2*ng)-ng+0.5)*self.dx + self.xl = xmin + (numpy.arange(nx+2*ng)-ng)*self.dx + self.xr = xmin + (numpy.arange(nx+2*ng)-ng+1.0)*self.dx + + # storage for the solution + # These are the modes of the solution at each point, so the + # first index is the mode + # NO! These are actually the *nodal* values + self.u = numpy.zeros((self.m+1, nx+2*ng), dtype=numpy.float64) + + # Need the Gauss-Lobatto nodes and weights in the reference element + GL = quadpy.line_segment.GaussLobatto(m+1) + self.nodes = GL.points + self.weights = GL.weights + # To go from modal to nodal we need the Vandermonde matrix + self.V = legendre.legvander(self.nodes, m) + c = numpy.eye(m+1) + # Orthonormalize + for p in range(m+1): + self.V[:, p] /= numpy.sqrt(2/(2*p+1)) + c[p, p] /= numpy.sqrt(2/(2*p+1)) + self.V_inv = numpy.linalg.inv(self.V) + self.M = numpy.linalg.inv(self.V @ self.V.T) + self.M_inv = self.V @ self.V.T + # Derivatives of Legendre polynomials lead to derivatives of V + dV = legendre.legval(self.nodes, + legendre.legder(c)).T + self.D = dV @ self.V_inv + # Stiffness matrix for the interior flux + self.S = self.M @ self.D + + # Nodes in the computational coordinates + self.all_nodes = numpy.zeros((m+1)*(nx+2*ng), dtype=numpy.float64) + self.all_nodes_per_node = numpy.zeros_like(self.u) + for i in range(nx+2*ng): + self.all_nodes[(m+1)*i:(m+1)*(i+1)] = (self.x[i] + + self.nodes * self.dx / 2) + self.all_nodes_per_node[:, i] = (self.x[i] + + self.nodes * self.dx / 2) + + def modal_to_nodal(self, modal): + for i in range(self.nx+2*self.ng): + self.u[:, i] = self.V @ modal[:, i] + return self.u + + def nodal_to_modal(self): + modal = numpy.zeros_like(self.u) + for i in range(self.nx+2*self.ng): + modal[:, i] = self.V_inv @ self.u[:, i] + return modal + + def plotting_data(self): + return (self.all_nodes, + self.u.ravel(order='F')) + + def plotting_data_high_order(self, npoints=50): + assert npoints > 2 + p_nodes = numpy.zeros(npoints*(self.nx+2*self.ng), dtype=numpy.float64) + p_data = numpy.zeros_like(p_nodes) + for i in range(self.nx+2*self.ng): + p_nodes[npoints*i:npoints*(i+1)] = numpy.linspace(self.xl[i], + self.xr[i], + npoints) + modal = self.V_inv @ self.u[:, i] + for p in range(self.m+1): + modal[p] /= numpy.sqrt(2/(2*p+1)) + scaled_x = 2 * (p_nodes[npoints*i:npoints*(i+1)] - + self.x[i]) / self.dx + p_data[npoints*i:npoints*(i+1)] = legendre.legval(scaled_x, modal) + return p_nodes, p_data + + def scratch_array(self): + """ return a scratch array dimensioned for our grid """ + return numpy.zeros((self.m+1, self.nx+2*self.ng), dtype=numpy.float64) + + def fill_BCs(self): + """ fill all single ghostcell with periodic boundary conditions """ + + for n in range(self.ng): + # left boundary + self.u[:, self.ilo-1-n] = self.u[:, self.ihi-n] + # right boundary + self.u[:, self.ihi+1+n] = self.u[:, self.ilo+n] + + def norm(self, e): + """ + Return the norm of quantity e which lives on the grid. + + This is the 'broken norm': the quantity is integrated over each + individual element using Gauss-Lobatto quadrature (as we have those + nodes and weights), and the 2-norm of the result is then returned. + """ + if not numpy.allclose(e.shape, self.all_nodes_per_node.shape): + return None + + # This is actually a pointwise norm, not quadrature'd + return numpy.sqrt(self.dx*numpy.sum(e[:, self.ilo:self.ihi+1]**2)) + + +@jit +def minmod3(a1, a2, a3): + """ + Utility function that does minmod on three inputs + """ + signs1 = a1 * a2 + signs2 = a1 * a3 + signs3 = a2 * a3 + same_sign = numpy.logical_and(numpy.logical_and(signs1 > 0, + signs2 > 0), + signs3 > 0) + minmod = numpy.min(numpy.abs(numpy.vstack((a1, a2, a3))), + axis=0) * numpy.sign(a1) + + return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) + + +@jit +def limit(g, limiter): + """ + After evolution, limit the slopes. + """ + + # Limiting! + + if limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[1:, :] = 0 + a_cell = (g.V @ a_zeromode)[0, :] + a_minus = g.u[0, :] + a_plus = g.u[-1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[i+1, 1:-1] + a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) + a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor + did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], + updated_mode[1:-1]) + a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + numpy.logical_not(did_it_limit)) + # Get back nodal values + limited_a = g.modal_to_nodal(a_modal) + + else: + limited_a = g.u + + return limited_a + + +class Simulation(object): + + def __init__(self, grid, C=0.8, limiter=None): + self.grid = grid + self.t = 0.0 # simulation time + self.C = C # CFL number + self.limiter = limiter # What it says. + + def init_cond(self, type="sine"): + """ initialize the data """ + if type == "tophat": + def init_u(x): + return numpy.where(numpy.logical_and(x >= 0.333, + x <= 0.666), + numpy.ones_like(x), + numpy.zeros_like(x)) + elif type == "sine": + def init_u(x): + return numpy.where(numpy.logical_and(x >= 0.333, + x <= 0.666), + numpy.ones_like(x) + + 0.5 * numpy.sin( + 2.0*numpy.pi* + (x-0.333)/0.333), + numpy.ones_like(x)) + elif type == "smooth_sine": + def init_u(x): + return numpy.sin(2.0 * numpy.pi * x / + (self.grid.xmax - self.grid.xmin)) + elif type == "gaussian": + def init_u(x): + local_xl = x - self.grid.dx/2 + local_xr = x + self.grid.dx/2 + al = 1.0 + numpy.exp(-60.0*(local_xl - 0.5)**2) + ar = 1.0 + numpy.exp(-60.0*(local_xr - 0.5)**2) + ac = 1.0 + numpy.exp(-60.0*(x - 0.5)**2) + + return (1./6.)*(al + 4*ac + ar) + + self.grid.u = init_u(self.grid.all_nodes_per_node) + + def timestep(self): + """ return the advective timestep """ + return self.C*self.grid.dx/numpy.max(numpy.abs( + self.grid.u[:, self.grid.ilo:self.grid.ihi+1])) + + + def states(self): + """ compute the left and right interface states """ + + # Evaluate the nodal values at the domain edges + g = self.grid + + # Extract nodal values + + al = numpy.zeros(g.nx+2*g.ng) + ar = numpy.zeros(g.nx+2*g.ng) + + # i is looping over interfaces, so al is the right edge of the left + # element, etc. + for i in range(g.ilo, g.ihi+2): + al[i] = g.u[-1, i-1] + ar[i] = g.u[0, i] + + return al, ar + + def limit(self): + """ + After evolution, limit the slopes. + """ + + # Evaluate the nodal values at the domain edges + g = self.grid + + # Limiting! + + if self.limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[1:, :] = 0 + a_cell = (g.V @ a_zeromode)[0, :] + a_minus = g.u[0, :] + a_plus = g.u[-1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[i+1, 1:-1] + a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) + a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor + did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], + updated_mode[1:-1]) + a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + numpy.logical_not(did_it_limit)) + # Get back nodal values + g.u = g.modal_to_nodal(a_modal) + + return g.u + + def riemann(self, ul, ur, alpha): + """ + Riemann problem for Burgers using Lax-Friedrichs + """ + + return ((ul**2 + ur**2)/2 - (ur - ul)*alpha)/2 + + def rk_substep(self, dt): + """ + Take a single RK substep + """ + g = self.grid + g.fill_BCs() + rhs = g.scratch_array() + + # Integrate flux over element + f = 0.5 * g.u**2 + interior_f = g.S.T @ f + # Use Riemann solver to get fluxes between elements + boundary_f = self.riemann(*self.states(), 2*numpy.max(numpy.abs(g.u))) + rhs = interior_f + rhs[0, 1:-1] += boundary_f[1:-1] + rhs[-1, 1:-1] -= boundary_f[2:] + + # Multiply by mass matrix (inverse). + rhs_i = 2 / g.dx * g.M_inv @ rhs + + return rhs_i + + def evolve(self, tmax): + """ evolve Burgers equation using RK3 (SSP) """ + self.t = 0.0 + g = self.grid + + # main evolution loop +# pbar = tqdm.tqdm(total=100) + while self.t < tmax: +# pbar.update(100*self.t/tmax) + # fill the boundary conditions + g.fill_BCs() + + # get the timestep + dt = self.timestep() + + if self.t + dt > tmax: + dt = tmax - self.t + + # RK3 - SSP + # Store the data at the start of the step + a_start = g.u.copy() + k1 = dt * self.rk_substep(dt) + g.u = a_start + k1 + g.u = self.limit() + g.fill_BCs() + a1 = g.u.copy() + k2 = dt * self.rk_substep(dt) + g.u = (3 * a_start + a1 + k2) / 4 + g.u = self.limit() + g.fill_BCs() + a2 = g.u.copy() + k3 = dt * self.rk_substep(dt) + g.u = (a_start + 2 * a2 + 2 * k3) / 3 + g.u = self.limit() + g.fill_BCs() + + self.t += dt +# pbar.close() + +if __name__ == "__main__": + + # Runs with limiter + nx = 128 + ng = 1 + xmin = 0 + xmax = 1 + ms = [1, 3, 5] + colors='brcy' + for i_m, m in enumerate(ms): + g = Grid1d(nx, ng, xmin, xmax, m) + s = Simulation(g, C=0.5/(2*m+1), limiter="moment") + s.init_cond("sine") + x, u = g.plotting_data() + x_start = x.copy() + u_start = u.copy() + s.evolve(0.2) + x_end, u_end = g.plotting_data() + pyplot.plot(x_end.copy(), u_end.copy(), f'{colors[i_m]}-', + label=rf'$m={m}$') + pyplot.plot(x_start, u_start, 'k--', alpha=0.5, label='Initial data') + pyplot.xlim(0, 1) + pyplot.legend() + pyplot.xlabel(r'$x$') + pyplot.ylabel(r'$u$') + pyplot.show() From fb9c0d7387b7fac1d9d085d0b3e62034448b1697 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 6 Mar 2019 21:13:25 +0000 Subject: [PATCH 29/33] Slight tidy up of Burgers --- burgers/dg_burgers.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/burgers/dg_burgers.py b/burgers/dg_burgers.py index e9a7bef..a9dea83 100644 --- a/burgers/dg_burgers.py +++ b/burgers/dg_burgers.py @@ -8,7 +8,7 @@ import matplotlib as mpl import quadpy from numba import jit -import tqdm +# import tqdm mpl.rcParams['mathtext.fontset'] = 'cm' mpl.rcParams['mathtext.rm'] = 'serif' @@ -133,6 +133,11 @@ def norm(self, e): return numpy.sqrt(self.dx*numpy.sum(e[:, self.ilo:self.ihi+1]**2)) +@jit +def burgers_flux(u): + return u**2/2 + + @jit def minmod3(a1, a2, a3): """ @@ -231,10 +236,10 @@ def init_u(x): def init_u(x): return numpy.where(numpy.logical_and(x >= 0.333, x <= 0.666), - numpy.ones_like(x) + + numpy.ones_like(x) + 0.5 * numpy.sin( - 2.0*numpy.pi* - (x-0.333)/0.333), + 2.0 * numpy.pi * + (x - 0.333) / 0.333), numpy.ones_like(x)) elif type == "smooth_sine": def init_u(x): @@ -257,7 +262,6 @@ def timestep(self): return self.C*self.grid.dx/numpy.max(numpy.abs( self.grid.u[:, self.grid.ilo:self.grid.ihi+1])) - def states(self): """ compute the left and right interface states """ @@ -341,7 +345,7 @@ def riemann(self, ul, ur, alpha): Riemann problem for Burgers using Lax-Friedrichs """ - return ((ul**2 + ur**2)/2 - (ur - ul)*alpha)/2 + return ((burgers_flux(ul) + burgers_flux(ur)) - (ur - ul)*alpha)/2 def rk_substep(self, dt): """ @@ -352,10 +356,10 @@ def rk_substep(self, dt): rhs = g.scratch_array() # Integrate flux over element - f = 0.5 * g.u**2 + f = burgers_flux(g.u) interior_f = g.S.T @ f # Use Riemann solver to get fluxes between elements - boundary_f = self.riemann(*self.states(), 2*numpy.max(numpy.abs(g.u))) + boundary_f = self.riemann(*self.states(), numpy.max(numpy.abs(g.u))) rhs = interior_f rhs[0, 1:-1] += boundary_f[1:-1] rhs[-1, 1:-1] -= boundary_f[2:] @@ -373,7 +377,7 @@ def evolve(self, tmax): # main evolution loop # pbar = tqdm.tqdm(total=100) while self.t < tmax: -# pbar.update(100*self.t/tmax) + # pbar.update(100*self.t/tmax) # fill the boundary conditions g.fill_BCs() @@ -404,6 +408,7 @@ def evolve(self, tmax): self.t += dt # pbar.close() + if __name__ == "__main__": # Runs with limiter @@ -412,7 +417,7 @@ def evolve(self, tmax): xmin = 0 xmax = 1 ms = [1, 3, 5] - colors='brcy' + colors = 'brcy' for i_m, m in enumerate(ms): g = Grid1d(nx, ng, xmin, xmax, m) s = Simulation(g, C=0.5/(2*m+1), limiter="moment") From 4e909f67dbe496485676073f9a55a70cbe941146 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 6 Mar 2019 21:13:42 +0000 Subject: [PATCH 30/33] Start DG Euler --- compressible/dg_euler.py | 493 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 493 insertions(+) create mode 100644 compressible/dg_euler.py diff --git a/compressible/dg_euler.py b/compressible/dg_euler.py new file mode 100644 index 0000000..ae1060b --- /dev/null +++ b/compressible/dg_euler.py @@ -0,0 +1,493 @@ +""" +Discontinuous Galerkin for the Euler equations. + +Just doing limiting on the conserved variables directly. +""" + +import numpy +from numpy.polynomial import legendre +from matplotlib import pyplot +import matplotlib as mpl +import quadpy +from numba import jit +# import tqdm +import riemann + +mpl.rcParams['mathtext.fontset'] = 'cm' +mpl.rcParams['mathtext.rm'] = 'serif' + + +class Grid1d(object): + + def __init__(self, nx, ng, xmin=0.0, xmax=1.0, m=3): + + assert m > 0 + + self.ncons = 3 # 1d hydro, 3 conserved variables + self.ng = ng + self.nx = nx + self.m = m + + self.xmin = xmin + self.xmax = xmax + + # 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 + (numpy.arange(nx+2*ng)-ng+0.5)*self.dx + self.xl = xmin + (numpy.arange(nx+2*ng)-ng)*self.dx + self.xr = xmin + (numpy.arange(nx+2*ng)-ng+1.0)*self.dx + + # storage for the solution + # These are the modes of the solution at each point, so the + # first index is the mode + # NO! These are actually the *nodal* values + self.u = numpy.zeros((self.ncons, self.m+1, nx+2*ng), + dtype=numpy.float64) + + # Need the Gauss-Lobatto nodes and weights in the reference element + GL = quadpy.line_segment.GaussLobatto(m+1) + self.nodes = GL.points + self.weights = GL.weights + # To go from modal to nodal we need the Vandermonde matrix + self.V = legendre.legvander(self.nodes, m) + c = numpy.eye(m+1) + # Orthonormalize + for p in range(m+1): + self.V[:, p] /= numpy.sqrt(2/(2*p+1)) + c[p, p] /= numpy.sqrt(2/(2*p+1)) + self.V_inv = numpy.linalg.inv(self.V) + self.M = numpy.linalg.inv(self.V @ self.V.T) + self.M_inv = self.V @ self.V.T + # Derivatives of Legendre polynomials lead to derivatives of V + dV = legendre.legval(self.nodes, + legendre.legder(c)).T + self.D = dV @ self.V_inv + # Stiffness matrix for the interior flux + self.S = self.M @ self.D + + # Nodes in the computational coordinates + self.all_nodes = numpy.zeros((m+1)*(nx+2*ng), dtype=numpy.float64) + self.all_nodes_per_node = numpy.zeros_like(self.u[0]) + for i in range(nx+2*ng): + self.all_nodes[(m+1)*i:(m+1)*(i+1)] = (self.x[i] + + self.nodes * self.dx / 2) + self.all_nodes_per_node[:, i] = (self.x[i] + + self.nodes * self.dx / 2) + + def modal_to_nodal(self, modal): + for n in range(self.ncons): + for i in range(self.nx+2*self.ng): + self.u[n, :, i] = self.V @ modal[n, :, i] + return self.u + + def nodal_to_modal(self): + modal = numpy.zeros_like(self.u) + for n in range(self.ncons): + for i in range(self.nx+2*self.ng): + modal[n, :, i] = self.V_inv @ self.u[n, :, i] + return modal + + def plotting_data(self): + u_plotting = numpy.zeros((self.ncons, (self.m+1)*(self.nx+2*self.ng))) + for n in range(self.ncons): + u_plotting[n, :] = self.u[n, :, :].ravel(order='F') + return (self.all_nodes, + u_plotting) + + def scratch_array(self): + """ return a scratch array dimensioned for our grid """ + return numpy.zeros((self.ncons, self.m+1, self.nx+2*self.ng), + dtype=numpy.float64) + + def fill_BCs(self): + """ fill all single ghostcell with outflow boundary conditions """ + + for n in range(self.ng): + # left boundary + self.u[:, :, self.ilo-1-n] = self.u[:, :, self.ilo] + # right boundary + self.u[:, :, self.ihi+1+n] = self.u[:, :, self.ihi] + + +@jit +def euler_flux(u, eos_gamma): + flux = numpy.zeros_like(u) + rho = u[0] + S = u[1] + E = u[2] + v = S / rho + p = (eos_gamma - 1) * (E - rho * v**2 / 2) + flux[0] = S + flux[1] = S * v + p + flux[2] = (E + p) * v + return flux + + +@jit +def minmod3(a1, a2, a3): + """ + Utility function that does minmod on three inputs + """ + signs1 = a1 * a2 + signs2 = a1 * a3 + signs3 = a2 * a3 + same_sign = numpy.logical_and(numpy.logical_and(signs1 > 0, + signs2 > 0), + signs3 > 0) + minmod = numpy.min(numpy.abs(numpy.vstack((a1, a2, a3))), + axis=0) * numpy.sign(a1) + + return numpy.where(same_sign, minmod, numpy.zeros_like(a1)) + + +@jit +def limit(g, limiter): + """ + After evolution, limit the slopes. + """ + + # Limiting! + + if limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[1:, :] = 0 + a_cell = (g.V @ a_zeromode)[0, :] + a_minus = g.u[0, :] + a_plus = g.u[-1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[i+1, 1:-1] + a2 = theta * (a_modal[i, 2:] - a_modal[i, 1:-1]) + a3 = theta * (a_modal[i, 1:-1] - a_modal[i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor + did_it_limit = numpy.isclose(a_modal[i+1, 1:-1], + updated_mode[1:-1]) + a_modal[i+1, limiting_todo] = updated_mode[limiting_todo] + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + numpy.logical_not(did_it_limit)) + # Get back nodal values + limited_a = g.modal_to_nodal(a_modal) + + else: + limited_a = g.u + + return limited_a + + +class Simulation(object): + + def __init__(self, grid, C=0.8, eos_gamma=1.4, limiter=None): + self.grid = grid + self.t = 0.0 # simulation time + self.C = C # CFL number + self.eos_gamma = eos_gamma + self.limiter = limiter # What it says. + + def init_cond(self, type="sod"): + """ initialize the data """ + if type == "sod": + rho_l = 1 + rho_r = 1 / 8 + v_l = 0 + v_r = 0 + p_l = 1 + p_r = 1 / 10 + S_l = rho_l * v_l + S_r = rho_r * v_r + e_l = p_l / rho_l / (self.eos_gamma - 1) + e_r = p_r / rho_r / (self.eos_gamma - 1) + E_l = rho_l * (e_l + v_l**2 / 2) + E_r = rho_r * (e_r + v_r**2 / 2) + x = self.grid.all_nodes_per_node + self.grid.u[0] = numpy.where(x < 0, + rho_l * numpy.ones_like(x), + rho_r * numpy.ones_like(x)) + self.grid.u[1] = numpy.where(x < 0, + S_l * numpy.ones_like(x), + S_r * numpy.ones_like(x)) + self.grid.u[2] = numpy.where(x < 0, + E_l * numpy.ones_like(x), + E_r * numpy.ones_like(x)) + elif type == "advection": + x = self.grid.all_nodes_per_node + rho_0 = 1e-3 + rho_1 = 1 + sigma = 0.1 + rho = rho_0 * numpy.ones_like(x) + rho += (rho_1 - rho_0) * numpy.exp(-(x-0.5)**2/sigma**2) + v = numpy.ones_like(x) + p = 1e-6 * numpy.ones_like(x) + S = rho * v + e = p / rho / (self.eos_gamma - 1) + E = rho * (e + v**2 / 2) + self.grid.u[0, :] = rho[:] + self.grid.u[1, :] = S[:] + self.grid.u[2, :] = E[:] + elif type == "double rarefaction": + rho_l = 1 + rho_r = 1 + v_l = -2 + v_r = 2 + p_l = 0.4 + p_r = 0.4 + S_l = rho_l * v_l + S_r = rho_r * v_r + e_l = p_l / rho_l / (self.eos_gamma - 1) + e_r = p_r / rho_r / (self.eos_gamma - 1) + E_l = rho_l * (e_l + v_l**2 / 2) + E_r = rho_r * (e_r + v_r**2 / 2) + x = self.grid.all_nodes_per_node + self.grid.u[0] = numpy.where(x < 0, + rho_l * numpy.ones_like(x), + rho_r * numpy.ones_like(x)) + self.grid.u[1] = numpy.where(x < 0, + S_l * numpy.ones_like(x), + S_r * numpy.ones_like(x)) + self.grid.u[2] = numpy.where(x < 0, + E_l * numpy.ones_like(x), + E_r * numpy.ones_like(x)) + + def max_lambda(self): + rho = self.grid.u[0] + v = self.grid.u[1] / rho + p = (self.eos_gamma - 1) * (self.grid.u[2, :] - rho * v**2 / 2) + cs = numpy.sqrt(self.eos_gamma * p / rho) + return numpy.max(numpy.abs(v) + cs) + + def timestep(self): + return self.C * self.grid.dx / self.max_lambda() + + def states(self): + """ compute the left and right interface states """ + + # Evaluate the nodal values at the domain edges + g = self.grid + + # Extract nodal values + al = numpy.zeros((g.ncons, g.nx+2*g.ng)) + ar = numpy.zeros((g.ncons, g.nx+2*g.ng)) + + # i is looping over interfaces, so al is the right edge of the left + # element, etc. + for i in range(g.ilo, g.ihi+2): + al[:, i] = g.u[:, -1, i-1] + ar[:, i] = g.u[:, 0, i] + + return al, ar + + def limit(self): + """ + After evolution, limit the slopes. + """ + + # Evaluate the nodal values at the domain edges + g = self.grid + + # Limiting! + + if self.limiter == "moment": + + # Limiting, using moment limiting (Hesthaven p 445-7) + theta = 2 + a_modal = g.nodal_to_modal() + # First, work out where limiting is needed + limiting_todo = numpy.ones(g.nx+2*g.ng, dtype=bool) + limiting_todo[:g.ng] = False + limiting_todo[-g.ng:] = False + # Get the cell average and the nodal values at the boundaries + a_zeromode = a_modal.copy() + a_zeromode[:, 1:, :] = 0 + for n in range(g.ncons): + a_cell = (g.V @ a_zeromode[n, :, :])[0, :] + a_minus = g.u[n, 0, :] + a_plus = g.u[n, -1, :] + # From the cell averages and boundary values we can construct + # alternate values at the boundaries + a_left = numpy.zeros(g.nx+2*g.ng) + a_right = numpy.zeros(g.nx+2*g.ng) + a_left[1:-1] = a_cell[1:-1] - minmod3(a_cell[1:-1] - a_minus[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + a_right[1:-1] = a_cell[1:-1] + minmod3(a_plus[1:-1] - a_cell[1:-1], + a_cell[1:-1] - a_cell[:-2], + a_cell[2:] - a_cell[1:-1]) + limiting_todo[1:-1] = numpy.logical_not( + numpy.logical_and(numpy.isclose(a_left[1:-1], + a_minus[1:-1]), + numpy.isclose(a_right[1:-1], + a_plus[1:-1]))) + # Now, do the limiting. Modify moments where needed, and as soon as + # limiting isn't needed, stop + updated_mode = numpy.zeros(g.nx+2*g.ng) + for i in range(g.m-1, -1, -1): + factor = numpy.sqrt((2*i+3)*(2*i+5)) + a1 = factor * a_modal[n, i+1, 1:-1] + a2 = theta * (a_modal[n, i, 2:] - a_modal[n, i, 1:-1]) + a3 = theta * (a_modal[n, i, 1:-1] - a_modal[n, i, :-2]) + updated_mode[1:-1] = minmod3(a1, a2, a3) / factor + did_it_limit = numpy.isclose(a_modal[n, i+1, 1:-1], + updated_mode[1:-1]) + a_modal[n, i+1, limiting_todo] = updated_mode[limiting_todo] + limiting_todo[1:-1] = numpy.logical_and(limiting_todo[1:-1], + numpy.logical_not(did_it_limit)) + # Get back nodal values + g.u = g.modal_to_nodal(a_modal) + + return g.u + + def riemann(self, ul, ur, alpha): + """ + Riemann problem for Burgers using Lax-Friedrichs + """ + fl = euler_flux(ul, self.eos_gamma) + fr = euler_flux(ur, self.eos_gamma) + return ((fl + fr) - (ur - ul)*alpha)/2 + + def rk_substep(self, dt): + """ + Take a single RK substep + """ + g = self.grid + g.fill_BCs() + rhs = g.scratch_array() + + # Integrate flux over element + f = euler_flux(g.u, self.eos_gamma) + interior_f = g.scratch_array() + for n in range(g.ncons): + interior_f[n, :, :] = g.S.T @ f[n, :, :] + # Use Riemann solver to get fluxes between elements + boundary_f = self.riemann(*self.states(), self.max_lambda()) + rhs = interior_f + rhs[:, 0, 1:-1] += boundary_f[:, 1:-1] + rhs[:, -1, 1:-1] -= boundary_f[:, 2:] + + # Multiply by mass matrix (inverse). + rhs_i = g.scratch_array() + for n in range(g.ncons): + rhs_i[n, :, :] = 2 / g.dx * g.M_inv @ rhs[n, :, :] + + return rhs_i + + def evolve(self, tmax): + """ evolve Burgers equation using RK3 (SSP) """ + self.t = 0.0 + g = self.grid + + # main evolution loop +# pbar = tqdm.tqdm(total=100) + while self.t < tmax: + # pbar.update(100*self.t/tmax) + # fill the boundary conditions + g.fill_BCs() + + # get the timestep + dt = self.timestep() + + if self.t + dt > tmax: + dt = tmax - self.t + + # RK3 - SSP + # Store the data at the start of the step + a_start = g.u.copy() + k1 = dt * self.rk_substep(dt) + g.u = a_start + k1 + g.u = self.limit() + g.fill_BCs() + a1 = g.u.copy() + k2 = dt * self.rk_substep(dt) + g.u = (3 * a_start + a1 + k2) / 4 + g.u = self.limit() + g.fill_BCs() + a2 = g.u.copy() + k3 = dt * self.rk_substep(dt) + g.u = (a_start + 2 * a2 + 2 * k3) / 3 + g.u = self.limit() + g.fill_BCs() + + self.t += dt +# pbar.close() + + +if __name__ == "__main__": + + # setup the problem -- Sod + left = riemann.State(p=1.0, u=0.0, rho=1.0) + right = riemann.State(p=0.1, u=0.0, rho=0.125) + + rp = riemann.RiemannProblem(left, right) + rp.find_star_state() + + x_e, rho_e, v_e, p_e = rp.sample_solution(0.2, 1024) + e_e = p_e / 0.4 / rho_e + x_e -= 0.5 + + # Runs with limiter + nx = 32 + ng = 1 + xmin = -0.5 + xmax = 0.5 + ms = [1, 3, 5] + colors = 'brcy' + for i_m, m in enumerate(ms): + g = Grid1d(nx, ng, xmin, xmax, m) + s = Simulation(g, C=0.5/(2*m+1), limiter="moment") + s.init_cond("sod") + s.evolve(0.2) + x, u = g.plotting_data() + rho = u[0, :] + v = u[1, :] / u[0, :] + e = (u[2, :] - rho * v**2 / 2) / rho + p = (s.eos_gamma - 1) * (u[2, :] - rho * v**2 / 2) + fig, axes = pyplot.subplots(2, 2) + axes[0, 0].plot(x[g.ilo:g.ihi+1], rho[g.ilo:g.ihi+1], 'bo') + axes[0, 0].plot(x_e, rho_e, 'k--') + axes[0, 1].plot(x[g.ilo:g.ihi+1], v[g.ilo:g.ihi+1], 'bo') + axes[0, 1].plot(x_e, v_e, 'k--') + axes[1, 0].plot(x[g.ilo:g.ihi+1], p[g.ilo:g.ihi+1], 'bo') + axes[1, 0].plot(x_e, p_e, 'k--') + axes[1, 1].plot(x[g.ilo:g.ihi+1], e[g.ilo:g.ihi+1], 'bo') + axes[1, 1].plot(x_e, e_e, 'k--') + axes[1, 0].set_xlabel(r"$x$") + axes[1, 1].set_xlabel(r"$x$") + axes[0, 0].set_ylabel(r"$\rho$") + axes[0, 1].set_ylabel(r"$v$") + axes[1, 0].set_ylabel(r"$p$") + axes[1, 1].set_ylabel(r"$e$") + for ax in axes.flatten(): + ax.set_xlim(-0.5, 0.5) + fig.tight_layout() + pyplot.show() + From 3307d99ff10f7f18da44b95e53035db0e80f3354 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Thu, 7 Mar 2019 17:53:19 +0000 Subject: [PATCH 31/33] DG for Euler --- compressible/dg_euler.py | 58 +++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 24 deletions(-) diff --git a/compressible/dg_euler.py b/compressible/dg_euler.py index ae1060b..17cb8ac 100644 --- a/compressible/dg_euler.py +++ b/compressible/dg_euler.py @@ -441,8 +441,7 @@ def evolve(self, tmax): # pbar.close() -if __name__ == "__main__": - +def plot_sod(nx, filename=None): # setup the problem -- Sod left = riemann.State(p=1.0, u=0.0, rho=1.0) right = riemann.State(p=0.1, u=0.0, rho=0.125) @@ -455,15 +454,19 @@ def evolve(self, tmax): x_e -= 0.5 # Runs with limiter - nx = 32 ng = 1 xmin = -0.5 xmax = 0.5 ms = [1, 3, 5] - colors = 'brcy' + colors = 'bry' + fig, axes = pyplot.subplots(2, 2) + axes[0, 0].plot(x_e, rho_e, 'k-', linewidth=2) + axes[0, 1].plot(x_e, v_e, 'k-', linewidth=2, label="Exact") + axes[1, 0].plot(x_e, p_e, 'k-', linewidth=2) + axes[1, 1].plot(x_e, e_e, 'k-', linewidth=2) for i_m, m in enumerate(ms): g = Grid1d(nx, ng, xmin, xmax, m) - s = Simulation(g, C=0.5/(2*m+1), limiter="moment") + s = Simulation(g, C=0.1/(2*m+1), limiter="moment") s.init_cond("sod") s.evolve(0.2) x, u = g.plotting_data() @@ -471,23 +474,30 @@ def evolve(self, tmax): v = u[1, :] / u[0, :] e = (u[2, :] - rho * v**2 / 2) / rho p = (s.eos_gamma - 1) * (u[2, :] - rho * v**2 / 2) - fig, axes = pyplot.subplots(2, 2) - axes[0, 0].plot(x[g.ilo:g.ihi+1], rho[g.ilo:g.ihi+1], 'bo') - axes[0, 0].plot(x_e, rho_e, 'k--') - axes[0, 1].plot(x[g.ilo:g.ihi+1], v[g.ilo:g.ihi+1], 'bo') - axes[0, 1].plot(x_e, v_e, 'k--') - axes[1, 0].plot(x[g.ilo:g.ihi+1], p[g.ilo:g.ihi+1], 'bo') - axes[1, 0].plot(x_e, p_e, 'k--') - axes[1, 1].plot(x[g.ilo:g.ihi+1], e[g.ilo:g.ihi+1], 'bo') - axes[1, 1].plot(x_e, e_e, 'k--') - axes[1, 0].set_xlabel(r"$x$") - axes[1, 1].set_xlabel(r"$x$") - axes[0, 0].set_ylabel(r"$\rho$") - axes[0, 1].set_ylabel(r"$v$") - axes[1, 0].set_ylabel(r"$p$") - axes[1, 1].set_ylabel(r"$e$") - for ax in axes.flatten(): - ax.set_xlim(-0.5, 0.5) - fig.tight_layout() - pyplot.show() + axes[0, 0].plot(x, rho, f'{colors[i_m]}--') + axes[0, 1].plot(x, v, f'{colors[i_m]}--', label=fr"$m={m}$") + axes[1, 0].plot(x, p, f'{colors[i_m]}--') + axes[1, 1].plot(x, e, f'{colors[i_m]}--') + axes[1, 0].set_xlabel(r"$x$") + axes[1, 1].set_xlabel(r"$x$") + axes[0, 0].set_ylabel(r"$\rho$") + axes[0, 1].set_ylabel(r"$v$") + axes[1, 0].set_ylabel(r"$p$") + axes[1, 1].set_ylabel(r"$e$") + axes[0, 0].set_title("Discontinuous Galerkin") + axes[0, 1].set_title(rf"$N={nx}$, Sod problem") + for ax in axes.flatten(): + ax.set_xlim(-0.5, 0.5) + fig.tight_layout() + lgd = axes[0, 1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) + if filename: + fig.savefig(filename, + bbox_extra_artists=(lgd,), bbox_inches='tight') + pyplot.show() + +if __name__ == "__main__": + + # Runs with limiter + for nx in [32, 128]: + plot_sod(nx, f'dg_sod_{nx}.pdf') From f896628ef6cd6dd0121417a87894966c1478a0e4 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Fri, 8 Mar 2019 11:28:50 +0000 Subject: [PATCH 32/33] Shock entropy interaction problem --- compressible/dg_euler.py | 87 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 81 insertions(+), 6 deletions(-) diff --git a/compressible/dg_euler.py b/compressible/dg_euler.py index 17cb8ac..41b7fbc 100644 --- a/compressible/dg_euler.py +++ b/compressible/dg_euler.py @@ -10,7 +10,7 @@ import matplotlib as mpl import quadpy from numba import jit -# import tqdm +import tqdm import riemann mpl.rcParams['mathtext.fontset'] = 'cm' @@ -279,6 +279,30 @@ def init_cond(self, type="sod"): E_l * numpy.ones_like(x), E_r * numpy.ones_like(x)) + elif type == "shock-entropy": + x = self.grid.all_nodes_per_node + rho_l = 3.857143 * numpy.ones_like(x) + rho_r = 0.2 * numpy.sin(numpy.pi*x) + 1 + v_l = 2.629369 * numpy.ones_like(x) + v_r = 0 * numpy.ones_like(x) + p_l = 10.33333 * numpy.ones_like(x) + p_r = 1 * numpy.ones_like(x) + S_l = rho_l * v_l + S_r = rho_r * v_r + e_l = p_l / rho_l / (self.eos_gamma - 1) + e_r = p_r / rho_r / (self.eos_gamma - 1) + E_l = rho_l * (e_l + v_l**2 / 2) + E_r = rho_r * (e_r + v_r**2 / 2) + self.grid.u[0] = numpy.where(x < -4, + rho_l, + rho_r) + self.grid.u[1] = numpy.where(x < -4, + S_l, + S_r) + self.grid.u[2] = numpy.where(x < -4, + E_l, + E_r) + def max_lambda(self): rho = self.grid.u[0] v = self.grid.u[1] / rho @@ -407,9 +431,15 @@ def evolve(self, tmax): g = self.grid # main evolution loop -# pbar = tqdm.tqdm(total=100) + # Estimate number of timesteps + nt_estimate = int(tmax / self.timestep() * 10) + prev_it = 0 + pbar = tqdm.tqdm(total=nt_estimate) while self.t < tmax: - # pbar.update(100*self.t/tmax) + curr_it = int(nt_estimate * self.t / tmax) + pbar.update(curr_it - prev_it) + prev_it = curr_it +# print(self.t) # fill the boundary conditions g.fill_BCs() @@ -438,7 +468,7 @@ def evolve(self, tmax): g.fill_BCs() self.t += dt -# pbar.close() + pbar.close() def plot_sod(nx, filename=None): @@ -496,8 +526,53 @@ def plot_sod(nx, filename=None): pyplot.show() +def plot_shock_entropy(nx, filename=None): + # Runs with limiter + ng = 1 + xmin = -5 + xmax = 5 + ms = [1, 3, 5] + colors = 'bry' + fig, axes = pyplot.subplots(2, 2) + for i_m, m in enumerate(ms): + g = Grid1d(nx, ng, xmin, xmax, m) + s = Simulation(g, C=0.1/(2*m+1), limiter="moment") + s.init_cond("shock-entropy") + print(f"Evolving m={m}") + s.evolve(1.8) + x, u = g.plotting_data() + rho = u[0, :] + v = u[1, :] / u[0, :] + e = (u[2, :] - rho * v**2 / 2) / rho + p = (s.eos_gamma - 1) * (u[2, :] - rho * v**2 / 2) + axes[0, 0].plot(x, rho, f'{colors[i_m]}--') + axes[0, 1].plot(x, v, f'{colors[i_m]}--', label=fr"$m={m}$") + axes[1, 0].plot(x, p, f'{colors[i_m]}--') + axes[1, 1].plot(x, e, f'{colors[i_m]}--') + axes[1, 0].set_xlabel(r"$x$") + axes[1, 1].set_xlabel(r"$x$") + axes[0, 0].set_ylabel(r"$\rho$") + axes[0, 1].set_ylabel(r"$v$") + axes[1, 0].set_ylabel(r"$p$") + axes[1, 1].set_ylabel(r"$e$") + axes[0, 0].set_title("Discontinuous Galerkin") + axes[0, 1].set_title(rf"$N={nx}$, Sod problem") + for ax in axes.flatten(): + ax.set_xlim(xmin, xmax) + fig.tight_layout() + lgd = axes[0, 1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) + if filename: + fig.savefig(filename, + bbox_extra_artists=(lgd,), bbox_inches='tight') + pyplot.show() + + if __name__ == "__main__": # Runs with limiter - for nx in [32, 128]: - plot_sod(nx, f'dg_sod_{nx}.pdf') +# for nx in [32, 128]: +# plot_sod(nx, f'dg_sod_{nx}.pdf') + + for nx in [128, 256]: + plot_shock_entropy(nx) + From 44e3a3bebe5269b0dd57ab780bb1cf4c73aab9b4 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Fri, 8 Mar 2019 15:24:58 +0000 Subject: [PATCH 33/33] More on shock entropy problem --- compressible/dg_euler.py | 72 ++++++++++++++++++++++++++------------ compressible/weno_euler.py | 33 ++++++++++++++++- 2 files changed, 82 insertions(+), 23 deletions(-) diff --git a/compressible/dg_euler.py b/compressible/dg_euler.py index 41b7fbc..8fa0f96 100644 --- a/compressible/dg_euler.py +++ b/compressible/dg_euler.py @@ -12,6 +12,7 @@ from numba import jit import tqdm import riemann +import weno_euler mpl.rcParams['mathtext.fontset'] = 'cm' mpl.rcParams['mathtext.rm'] = 'serif' @@ -278,7 +279,6 @@ def init_cond(self, type="sod"): self.grid.u[2] = numpy.where(x < 0, E_l * numpy.ones_like(x), E_r * numpy.ones_like(x)) - elif type == "shock-entropy": x = self.grid.all_nodes_per_node rho_l = 3.857143 * numpy.ones_like(x) @@ -526,29 +526,52 @@ def plot_sod(nx, filename=None): pyplot.show() -def plot_shock_entropy(nx, filename=None): +def weno_reference_shock_entropy(nx=512): + """ + A reference solution for the shock-entropy problem. + + Note that this is slow: 20 minutes or so using the standard parameters, + and obviously more resolution would be better. + """ + order = 4 + ng = order+2 + C = 0.5 + xmin = -5 + xmax = 5 + g = weno_euler.Grid1d(nx, ng, xmin, xmax, bc="outflow") + s = weno_euler.WENOSimulation(g, C, order) + s.init_cond("shock-entropy") + s.evolve(1.8, reconstruction='characteristic') + return g.x, g.q + +def plot_shock_entropy(nx, m, x_ref, u_ref, filename=None): # Runs with limiter ng = 1 xmin = -5 xmax = 5 - ms = [1, 3, 5] - colors = 'bry' fig, axes = pyplot.subplots(2, 2) - for i_m, m in enumerate(ms): - g = Grid1d(nx, ng, xmin, xmax, m) - s = Simulation(g, C=0.1/(2*m+1), limiter="moment") - s.init_cond("shock-entropy") - print(f"Evolving m={m}") - s.evolve(1.8) - x, u = g.plotting_data() - rho = u[0, :] - v = u[1, :] / u[0, :] - e = (u[2, :] - rho * v**2 / 2) / rho - p = (s.eos_gamma - 1) * (u[2, :] - rho * v**2 / 2) - axes[0, 0].plot(x, rho, f'{colors[i_m]}--') - axes[0, 1].plot(x, v, f'{colors[i_m]}--', label=fr"$m={m}$") - axes[1, 0].plot(x, p, f'{colors[i_m]}--') - axes[1, 1].plot(x, e, f'{colors[i_m]}--') + rho_ref = u_ref[0, :] + v_ref = u_ref[1, :] / u_ref[0, :] + e_ref = (u_ref[2, :] - rho_ref * v_ref**2 / 2) / rho_ref + p_ref = 0.4 * (u_ref[2, :] - rho_ref * v_ref**2 / 2) + axes[0, 0].plot(x_ref, rho_ref, 'k--') + axes[0, 1].plot(x_ref, v_ref, 'k--', label="Reference") + axes[1, 0].plot(x_ref, p_ref, 'k--') + axes[1, 1].plot(x_ref, e_ref, 'k--') + # Evolve DG + g = Grid1d(nx, ng, xmin, xmax, m) + s = Simulation(g, C=0.1/(2*m+1), limiter="moment") + s.init_cond("shock-entropy") + s.evolve(1.8) + x, u = g.plotting_data() + rho = u[0, :] + v = u[1, :] / u[0, :] + e = (u[2, :] - rho * v**2 / 2) / rho + p = (s.eos_gamma - 1) * (u[2, :] - rho * v**2 / 2) + axes[0, 0].plot(x, rho, 'b-') + axes[0, 1].plot(x, v, 'b-', label=fr"$m={m}$") + axes[1, 0].plot(x, p, 'b-') + axes[1, 1].plot(x, e, 'b-') axes[1, 0].set_xlabel(r"$x$") axes[1, 1].set_xlabel(r"$x$") axes[0, 0].set_ylabel(r"$\rho$") @@ -556,7 +579,7 @@ def plot_shock_entropy(nx, filename=None): axes[1, 0].set_ylabel(r"$p$") axes[1, 1].set_ylabel(r"$e$") axes[0, 0].set_title("Discontinuous Galerkin") - axes[0, 1].set_title(rf"$N={nx}$, Sod problem") + axes[0, 1].set_title(rf"$N={nx}$, shock-entropy problem") for ax in axes.flatten(): ax.set_xlim(xmin, xmax) fig.tight_layout() @@ -573,6 +596,11 @@ def plot_shock_entropy(nx, filename=None): # for nx in [32, 128]: # plot_sod(nx, f'dg_sod_{nx}.pdf') - for nx in [128, 256]: - plot_shock_entropy(nx) + # Reference solution for shock-entropy using WENO + x_ref, u_ref = weno_reference_shock_entropy() + + for m in [1, 3, 5]: + for nx in [128, 256]: + plot_shock_entropy(nx, m, x_ref, u_ref, + f'dg_shock_entropy_m{m}_N{nx}.pdf') diff --git a/compressible/weno_euler.py b/compressible/weno_euler.py index fd7575d..c33307e 100644 --- a/compressible/weno_euler.py +++ b/compressible/weno_euler.py @@ -3,6 +3,7 @@ from matplotlib import pyplot import weno_coefficients import riemann +import tqdm class Grid1d(object): @@ -172,6 +173,29 @@ def init_cond(self, type="sod"): self.grid.q[2] = numpy.where(self.grid.x < 0, E_l * numpy.ones_like(self.grid.x), E_r * numpy.ones_like(self.grid.x)) + elif type == "shock-entropy": + x = self.grid.x + rho_l = 3.857143 * numpy.ones_like(x) + rho_r = 0.2 * numpy.sin(numpy.pi*x) + 1 + v_l = 2.629369 * numpy.ones_like(x) + v_r = 0 * numpy.ones_like(x) + p_l = 10.33333 * numpy.ones_like(x) + p_r = 1 * numpy.ones_like(x) + S_l = rho_l * v_l + S_r = rho_r * v_r + e_l = p_l / rho_l / (self.eos_gamma - 1) + e_r = p_r / rho_r / (self.eos_gamma - 1) + E_l = rho_l * (e_l + v_l**2 / 2) + E_r = rho_r * (e_r + v_r**2 / 2) + self.grid.q[0] = numpy.where(x < -4, + rho_l, + rho_r) + self.grid.q[1] = numpy.where(x < -4, + S_l, + S_r) + self.grid.q[2] = numpy.where(x < -4, + E_l, + E_r) def max_lambda(self): @@ -297,7 +321,14 @@ def evolve(self, tmax, reconstruction = 'componentwise'): stepper = self.rk_substep_characteristic # main evolution loop + # Estimate number of timesteps + nt_estimate = int(tmax / self.timestep() * 10) + prev_it = 0 + pbar = tqdm.tqdm(total=nt_estimate) while self.t < tmax: + curr_it = int(nt_estimate * self.t / tmax) + pbar.update(curr_it - prev_it) + prev_it = curr_it # fill the boundary conditions g.fill_BCs() @@ -331,7 +362,7 @@ def evolve(self, tmax, reconstruction = 'componentwise'): self.t += dt # print("t=", self.t) - + pbar.close() if __name__ == "__main__":