From c26de495738abe1d4d67a3e2e657f69eb775cdd6 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Fri, 26 Jul 2019 15:42:10 -0400 Subject: [PATCH 01/47] here --- .../.acoustics_2d_inclusions.py.swp | Bin 0 -> 1024 bytes .../acoustics_2d_mapped/inclusion_mapping.py | 297 ++++++++++++++++++ examples/acoustics_2d_mapped/setplot.py | 297 ++++++++++++++++++ examples/shallow_2d/plot.py | 80 +++++ src/pyclaw/classic/solver.py | 167 +++++----- 5 files changed, 763 insertions(+), 78 deletions(-) create mode 100644 examples/acoustics_2d_mapped/.acoustics_2d_inclusions.py.swp create mode 100755 examples/acoustics_2d_mapped/inclusion_mapping.py create mode 100755 examples/acoustics_2d_mapped/setplot.py create mode 100644 examples/shallow_2d/plot.py diff --git a/examples/acoustics_2d_mapped/.acoustics_2d_inclusions.py.swp b/examples/acoustics_2d_mapped/.acoustics_2d_inclusions.py.swp new file mode 100644 index 0000000000000000000000000000000000000000..e721fc9b88fedc954a03e52adc18333d7ee55f5c GIT binary patch literal 1024 zcmYc?$V<%2S1{7E)H7y40#X_b3|XZqi5W;@xR_vdx*-{vdD#JpDY_9x#s*lGCMM^X n7MElu7snf=#AoIu=ad#_=I0gb6;xtXJ}NjG0;3^7;}8G<4cHcq literal 0 HcmV?d00001 diff --git a/examples/acoustics_2d_mapped/inclusion_mapping.py b/examples/acoustics_2d_mapped/inclusion_mapping.py new file mode 100755 index 000000000..8981787a5 --- /dev/null +++ b/examples/acoustics_2d_mapped/inclusion_mapping.py @@ -0,0 +1,297 @@ +#!/usr/bin/env python +# encoding: utf-8 +r""" +Two-dimensional variable-coefficient acoustics on a mapped grid +=============================================================== + +Solve the variable-coefficient acoustics equations in 2D: + +.. math:: + p_t + K(x,y) (u_x + v_y) & = 0 \\ + u_t + p_x / \rho(x,y) & = 0 \\ + v_t + p_y / \rho(x,y) & = 0. + +Here p is the pressure, (u,v) is the velocity, :math:`K(x,y)` is the bulk modulus, +and :math:`\rho(x,y)` is the density. + +This example shows how to solve a problem with variable coefficients on a mapped grid. +The domain contains circular inclusions with different acoustic properties. +""" +from __future__ import absolute_import +import numpy as np +from six.moves import range + +# Circle radius, square radius, circle center: +# ((r1, r2), (x0, y0)) +circles = ( ((0.15,0.204),( 0.45,0.3)), + ((0.15,0.204),( 0.7, 0.75))) + +impedance = (12.,10.) +sound_speed = (0.3,1.5) + +def inclusion_mapping(xc, yc): + """Apply mapping from square to circle for each circle specified. + Leave rest of grid cartesian. + """ + xp = xc+0. + yp = yc+0. + + for circle in circles: + circle_center = circle[1] + r1 = circle[0][0] + r2 = circle[0][1] + + x0, y0 = circle_center + xdist = np.abs(xc-x0) + ydist = np.abs(yc-y0) + insquare = np.where(np.maximum(xdist,ydist)<=r2) # Square section of grid to be deformed + xc0 = (xc[insquare]-x0)/r2 + yc0 = (yc[insquare]-y0)/r2 + + xc1 = np.abs(xc0) + yc1 = np.abs(yc0) + d = np.maximum(xc1, yc1) + d = np.maximum(d,1.e-10) + d = np.minimum(d, 0.99999) + d1 = d*r2/np.sqrt(2) + + # Pick mapping + #R = np.sqrt(2) * d1 + R = r1*np.ones(d1.shape) + R = r1**2 / (r2*d) + + # Modify d1 and R ouside circle to morph back to square: + ij = np.where(d>(r1/r2)) + d1[ij] = r1/np.sqrt(2) + (d[ij]-r1/r2)*(r2-r1/np.sqrt(2))/(1.-r1/r2) + R[ij] = r1 * ((1.-r1/r2) / (1.-d[ij]))**(r2/r1 + 0.5) + + xp2 = d1/d * xc1 + yp2 = d1/d * yc1 + center = d1 - np.sqrt(R**2 - d1**2) + + ij = np.where(xc1>=yc1) + xp2[ij] = center[ij] + np.sqrt(R[ij]**2 - yp2[ij]**2) + + ij = np.where(yc1>=xc1) + yp2[ij] = center[ij] + np.sqrt(R[ij]**2 - xp2[ij]**2) + + xp2 = np.sign(xc0) * xp2 + yp2 = np.sign(yc0) * yp2 + + xp[insquare] = x0 + xp2 + yp[insquare] = y0 + yp2 + + return xp, yp + +def compute_geometry(grid): + r"""Computes + a_x + a_y + length_ratio_left + b_x + b_y + length_ratio_bottom + cell_area + """ + + dx, dy = grid.delta + area_min = 1.e6 + area_max = 0.0 + + x_corners, y_corners = grid.p_nodes + + lower_left_y, lower_left_x = y_corners[:-1,:-1], x_corners[:-1,:-1] + upper_left_y, upper_left_x = y_corners[:-1,1: ], x_corners[:-1,1: ] + lower_right_y, lower_right_x = y_corners[1:,:-1], x_corners[1:,:-1] + upper_right_y, upper_right_x = y_corners[1:,1: ], x_corners[1:,1: ] + + a_x = upper_left_y - lower_left_y #upper left and lower left + a_y = -(upper_left_x - lower_left_x) + anorm = np.sqrt(a_x**2 + a_y**2) + a_x, a_y = a_x/anorm, a_y/anorm + length_ratio_left = anorm/dy + + b_x = -(lower_right_y - lower_left_y) #lower right and lower left + b_y = lower_right_x - lower_left_x + bnorm = np.sqrt(b_x**2 + b_y**2) + b_x, b_y = b_x/bnorm, b_y/bnorm + length_ratio_bottom = bnorm/dx + + area = 0*grid.c_centers[0] + area += 0.5 * (lower_left_y+upper_left_y)*(upper_left_x-lower_left_x) + area += 0.5 * (upper_left_y+upper_right_y)*(upper_right_x-upper_left_x) + area += 0.5 * (upper_right_y+lower_right_y)*(lower_right_x-upper_right_x) + area += 0.5 * (lower_right_y+lower_left_y)*(lower_left_x-lower_right_x) + area = area/(dx*dy) + area_min = min(area_min, np.min(area)) + area_max = max(area_max, np.max(area)) + + return a_x, a_y, length_ratio_left, b_x, b_y, length_ratio_bottom, area + +def incoming_square_wave(state,dim,t,qbc,auxbc,num_ghost): + """ + Incoming square wave at left boundary. + """ + if t<0.05: + s = 1.0 + else: + s = 0.0 + + for i in range(num_ghost): + reflect_ind = 2*num_ghost - i - 1 + alpha = auxbc[0,i,:] + beta = auxbc[1,i,:] + u_normal = alpha*qbc[1,reflect_ind,:] + beta*qbc[2,reflect_ind,:] + u_tangential = -beta*qbc[1,reflect_ind,:] + alpha*qbc[2,reflect_ind,:] + u_normal = 2.0*s - u_normal + qbc[0,i,:] = qbc[0,reflect_ind,:] + qbc[1,i,:] = alpha*u_normal - beta*u_tangential + qbc[2,i,:] = beta*u_normal + alpha*u_tangential + + +def setup(kernel_language='Fortran', use_petsc=False, outdir='./_output', + solver_type='classic', time_integrator='SSP104', lim_type=2, + num_output_times=20, disable_output=False, num_cells=200): + from clawpack import riemann + + if use_petsc: + import clawpack.petclaw as pyclaw + else: + from clawpack import pyclaw + + riemann_solver = riemann.acoustics_mapped_2D + + if solver_type=='classic': + solver=pyclaw.ClawSolver2D(riemann_solver) + solver.dimensional_split=False + solver.limiters = pyclaw.limiters.tvd.MC + elif solver_type=='sharpclaw': + solver=pyclaw.SharpClawSolver2D(riemann_solver) + solver.time_integrator=time_integrator + + solver.bc_lower[0]=pyclaw.BC.custom + solver.user_bc_lower = incoming_square_wave + solver.bc_upper[0]=pyclaw.BC.extrap + solver.bc_lower[1]=pyclaw.BC.extrap + solver.bc_upper[1]=pyclaw.BC.extrap + + solver.aux_bc_lower[0]=pyclaw.BC.extrap + solver.aux_bc_upper[0]=pyclaw.BC.extrap + solver.aux_bc_lower[1]=pyclaw.BC.extrap + solver.aux_bc_upper[1]=pyclaw.BC.extrap + + x = pyclaw.Dimension(0.,1.0,num_cells,name='x') + y = pyclaw.Dimension(0.,1.0,num_cells,name='y') + domain = pyclaw.Domain([x,y]) + + num_eqn = 3 + num_aux = 9 # geometry (7), impedance, sound speed + state = pyclaw.State(domain,num_eqn,num_aux) + state.grid.mapc2p = inclusion_mapping + + a_x, a_y, length_left, b_x, b_y, length_bottom, area = compute_geometry(state.grid) + + state.aux[0,:,:] = a_x + state.aux[1,:,:] = a_y + state.aux[2,:,:] = length_left + state.aux[3,:,:] = b_x + state.aux[4,:,:] = b_y + state.aux[5,:,:] = length_bottom + state.aux[6,:,:] = area + state.index_capa = 6 # aux[6,:,:] holds the capacity function + + grid = state.grid + xp, yp = grid.p_centers + state.aux[7,:,:] = 1.0 # Impedance + state.aux[8,:,:] = 1.0 # Sound speed + + for i, circle in enumerate(circles): + # Set impedance and sound speed in each inclusion + radius = circle[0][0] + x0, y0 = circle[1] + distance = np.sqrt( (xp-x0)**2 + (yp-y0)**2 ) + in_circle = np.where(distance <= radius) + state.aux[7][in_circle] = impedance[i] + state.aux[8][in_circle] = sound_speed[i] + + # Set initial condition + state.q[0,:,:] = 0. + state.q[1,:,:] = 0. + state.q[2,:,:] = 0. + + claw = pyclaw.Controller() + claw.keep_copy = True + if disable_output: + claw.output_format = None + claw.solution = pyclaw.Solution(state,domain) + claw.solver = solver + claw.outdir = outdir + claw.tfinal = 0.9 + claw.num_output_times = num_output_times + claw.write_aux_init = True + claw.setplot = setplot + if use_petsc: + claw.output_options = {'format':'binary'} + + return claw + + +def setplot(plotdata): + """ + Plot solution using VisClaw. + + This example shows how to mark an internal boundary on a 2D plot. + """ + + from clawpack.visclaw import colormaps + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.mapc2p = inclusion_mapping + + # Figure for pressure + plotfigure = plotdata.new_plotfigure(name='Pressure', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.title = 'Pressure' + plotaxes.scaled = True # so aspect ratio is 1 + plotaxes.afteraxes = plot_circles + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.contour_nlevels = 100 + plotitem.contour_min = -2.51 + plotitem.contour_max = 2.51 + plotitem.plot_var = 0 + + # Figure for x-velocity plot + plotfigure = plotdata.new_plotfigure(name='x-Velocity', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.title = 'u' + plotaxes.afteraxes = plot_circles + + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = 1 + plotitem.pcolor_cmap = 'RdBu' + plotitem.add_colorbar = True + plotitem.pcolor_cmin = -1.0 + plotitem.pcolor_cmax= 1.0 + + return plotdata + +def plot_circles(current_data): + import matplotlib.pyplot as plt + ax = plt.gca() + for circle in circles: + x0, y0 = circle[1] + radius = circle[0][0] + circle1=plt.Circle((x0,y0),radius,color='k',fill=False,lw=3) + ax.add_artist(circle1) + + +if __name__=="__main__": + import sys + from clawpack.pyclaw.util import run_app_from_main + output = run_app_from_main(setup,setplot) diff --git a/examples/acoustics_2d_mapped/setplot.py b/examples/acoustics_2d_mapped/setplot.py new file mode 100755 index 000000000..8981787a5 --- /dev/null +++ b/examples/acoustics_2d_mapped/setplot.py @@ -0,0 +1,297 @@ +#!/usr/bin/env python +# encoding: utf-8 +r""" +Two-dimensional variable-coefficient acoustics on a mapped grid +=============================================================== + +Solve the variable-coefficient acoustics equations in 2D: + +.. math:: + p_t + K(x,y) (u_x + v_y) & = 0 \\ + u_t + p_x / \rho(x,y) & = 0 \\ + v_t + p_y / \rho(x,y) & = 0. + +Here p is the pressure, (u,v) is the velocity, :math:`K(x,y)` is the bulk modulus, +and :math:`\rho(x,y)` is the density. + +This example shows how to solve a problem with variable coefficients on a mapped grid. +The domain contains circular inclusions with different acoustic properties. +""" +from __future__ import absolute_import +import numpy as np +from six.moves import range + +# Circle radius, square radius, circle center: +# ((r1, r2), (x0, y0)) +circles = ( ((0.15,0.204),( 0.45,0.3)), + ((0.15,0.204),( 0.7, 0.75))) + +impedance = (12.,10.) +sound_speed = (0.3,1.5) + +def inclusion_mapping(xc, yc): + """Apply mapping from square to circle for each circle specified. + Leave rest of grid cartesian. + """ + xp = xc+0. + yp = yc+0. + + for circle in circles: + circle_center = circle[1] + r1 = circle[0][0] + r2 = circle[0][1] + + x0, y0 = circle_center + xdist = np.abs(xc-x0) + ydist = np.abs(yc-y0) + insquare = np.where(np.maximum(xdist,ydist)<=r2) # Square section of grid to be deformed + xc0 = (xc[insquare]-x0)/r2 + yc0 = (yc[insquare]-y0)/r2 + + xc1 = np.abs(xc0) + yc1 = np.abs(yc0) + d = np.maximum(xc1, yc1) + d = np.maximum(d,1.e-10) + d = np.minimum(d, 0.99999) + d1 = d*r2/np.sqrt(2) + + # Pick mapping + #R = np.sqrt(2) * d1 + R = r1*np.ones(d1.shape) + R = r1**2 / (r2*d) + + # Modify d1 and R ouside circle to morph back to square: + ij = np.where(d>(r1/r2)) + d1[ij] = r1/np.sqrt(2) + (d[ij]-r1/r2)*(r2-r1/np.sqrt(2))/(1.-r1/r2) + R[ij] = r1 * ((1.-r1/r2) / (1.-d[ij]))**(r2/r1 + 0.5) + + xp2 = d1/d * xc1 + yp2 = d1/d * yc1 + center = d1 - np.sqrt(R**2 - d1**2) + + ij = np.where(xc1>=yc1) + xp2[ij] = center[ij] + np.sqrt(R[ij]**2 - yp2[ij]**2) + + ij = np.where(yc1>=xc1) + yp2[ij] = center[ij] + np.sqrt(R[ij]**2 - xp2[ij]**2) + + xp2 = np.sign(xc0) * xp2 + yp2 = np.sign(yc0) * yp2 + + xp[insquare] = x0 + xp2 + yp[insquare] = y0 + yp2 + + return xp, yp + +def compute_geometry(grid): + r"""Computes + a_x + a_y + length_ratio_left + b_x + b_y + length_ratio_bottom + cell_area + """ + + dx, dy = grid.delta + area_min = 1.e6 + area_max = 0.0 + + x_corners, y_corners = grid.p_nodes + + lower_left_y, lower_left_x = y_corners[:-1,:-1], x_corners[:-1,:-1] + upper_left_y, upper_left_x = y_corners[:-1,1: ], x_corners[:-1,1: ] + lower_right_y, lower_right_x = y_corners[1:,:-1], x_corners[1:,:-1] + upper_right_y, upper_right_x = y_corners[1:,1: ], x_corners[1:,1: ] + + a_x = upper_left_y - lower_left_y #upper left and lower left + a_y = -(upper_left_x - lower_left_x) + anorm = np.sqrt(a_x**2 + a_y**2) + a_x, a_y = a_x/anorm, a_y/anorm + length_ratio_left = anorm/dy + + b_x = -(lower_right_y - lower_left_y) #lower right and lower left + b_y = lower_right_x - lower_left_x + bnorm = np.sqrt(b_x**2 + b_y**2) + b_x, b_y = b_x/bnorm, b_y/bnorm + length_ratio_bottom = bnorm/dx + + area = 0*grid.c_centers[0] + area += 0.5 * (lower_left_y+upper_left_y)*(upper_left_x-lower_left_x) + area += 0.5 * (upper_left_y+upper_right_y)*(upper_right_x-upper_left_x) + area += 0.5 * (upper_right_y+lower_right_y)*(lower_right_x-upper_right_x) + area += 0.5 * (lower_right_y+lower_left_y)*(lower_left_x-lower_right_x) + area = area/(dx*dy) + area_min = min(area_min, np.min(area)) + area_max = max(area_max, np.max(area)) + + return a_x, a_y, length_ratio_left, b_x, b_y, length_ratio_bottom, area + +def incoming_square_wave(state,dim,t,qbc,auxbc,num_ghost): + """ + Incoming square wave at left boundary. + """ + if t<0.05: + s = 1.0 + else: + s = 0.0 + + for i in range(num_ghost): + reflect_ind = 2*num_ghost - i - 1 + alpha = auxbc[0,i,:] + beta = auxbc[1,i,:] + u_normal = alpha*qbc[1,reflect_ind,:] + beta*qbc[2,reflect_ind,:] + u_tangential = -beta*qbc[1,reflect_ind,:] + alpha*qbc[2,reflect_ind,:] + u_normal = 2.0*s - u_normal + qbc[0,i,:] = qbc[0,reflect_ind,:] + qbc[1,i,:] = alpha*u_normal - beta*u_tangential + qbc[2,i,:] = beta*u_normal + alpha*u_tangential + + +def setup(kernel_language='Fortran', use_petsc=False, outdir='./_output', + solver_type='classic', time_integrator='SSP104', lim_type=2, + num_output_times=20, disable_output=False, num_cells=200): + from clawpack import riemann + + if use_petsc: + import clawpack.petclaw as pyclaw + else: + from clawpack import pyclaw + + riemann_solver = riemann.acoustics_mapped_2D + + if solver_type=='classic': + solver=pyclaw.ClawSolver2D(riemann_solver) + solver.dimensional_split=False + solver.limiters = pyclaw.limiters.tvd.MC + elif solver_type=='sharpclaw': + solver=pyclaw.SharpClawSolver2D(riemann_solver) + solver.time_integrator=time_integrator + + solver.bc_lower[0]=pyclaw.BC.custom + solver.user_bc_lower = incoming_square_wave + solver.bc_upper[0]=pyclaw.BC.extrap + solver.bc_lower[1]=pyclaw.BC.extrap + solver.bc_upper[1]=pyclaw.BC.extrap + + solver.aux_bc_lower[0]=pyclaw.BC.extrap + solver.aux_bc_upper[0]=pyclaw.BC.extrap + solver.aux_bc_lower[1]=pyclaw.BC.extrap + solver.aux_bc_upper[1]=pyclaw.BC.extrap + + x = pyclaw.Dimension(0.,1.0,num_cells,name='x') + y = pyclaw.Dimension(0.,1.0,num_cells,name='y') + domain = pyclaw.Domain([x,y]) + + num_eqn = 3 + num_aux = 9 # geometry (7), impedance, sound speed + state = pyclaw.State(domain,num_eqn,num_aux) + state.grid.mapc2p = inclusion_mapping + + a_x, a_y, length_left, b_x, b_y, length_bottom, area = compute_geometry(state.grid) + + state.aux[0,:,:] = a_x + state.aux[1,:,:] = a_y + state.aux[2,:,:] = length_left + state.aux[3,:,:] = b_x + state.aux[4,:,:] = b_y + state.aux[5,:,:] = length_bottom + state.aux[6,:,:] = area + state.index_capa = 6 # aux[6,:,:] holds the capacity function + + grid = state.grid + xp, yp = grid.p_centers + state.aux[7,:,:] = 1.0 # Impedance + state.aux[8,:,:] = 1.0 # Sound speed + + for i, circle in enumerate(circles): + # Set impedance and sound speed in each inclusion + radius = circle[0][0] + x0, y0 = circle[1] + distance = np.sqrt( (xp-x0)**2 + (yp-y0)**2 ) + in_circle = np.where(distance <= radius) + state.aux[7][in_circle] = impedance[i] + state.aux[8][in_circle] = sound_speed[i] + + # Set initial condition + state.q[0,:,:] = 0. + state.q[1,:,:] = 0. + state.q[2,:,:] = 0. + + claw = pyclaw.Controller() + claw.keep_copy = True + if disable_output: + claw.output_format = None + claw.solution = pyclaw.Solution(state,domain) + claw.solver = solver + claw.outdir = outdir + claw.tfinal = 0.9 + claw.num_output_times = num_output_times + claw.write_aux_init = True + claw.setplot = setplot + if use_petsc: + claw.output_options = {'format':'binary'} + + return claw + + +def setplot(plotdata): + """ + Plot solution using VisClaw. + + This example shows how to mark an internal boundary on a 2D plot. + """ + + from clawpack.visclaw import colormaps + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.mapc2p = inclusion_mapping + + # Figure for pressure + plotfigure = plotdata.new_plotfigure(name='Pressure', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.title = 'Pressure' + plotaxes.scaled = True # so aspect ratio is 1 + plotaxes.afteraxes = plot_circles + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.contour_nlevels = 100 + plotitem.contour_min = -2.51 + plotitem.contour_max = 2.51 + plotitem.plot_var = 0 + + # Figure for x-velocity plot + plotfigure = plotdata.new_plotfigure(name='x-Velocity', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.title = 'u' + plotaxes.afteraxes = plot_circles + + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = 1 + plotitem.pcolor_cmap = 'RdBu' + plotitem.add_colorbar = True + plotitem.pcolor_cmin = -1.0 + plotitem.pcolor_cmax= 1.0 + + return plotdata + +def plot_circles(current_data): + import matplotlib.pyplot as plt + ax = plt.gca() + for circle in circles: + x0, y0 = circle[1] + radius = circle[0][0] + circle1=plt.Circle((x0,y0),radius,color='k',fill=False,lw=3) + ax.add_artist(circle1) + + +if __name__=="__main__": + import sys + from clawpack.pyclaw.util import run_app_from_main + output = run_app_from_main(setup,setplot) diff --git a/examples/shallow_2d/plot.py b/examples/shallow_2d/plot.py new file mode 100644 index 000000000..f46427875 --- /dev/null +++ b/examples/shallow_2d/plot.py @@ -0,0 +1,80 @@ +r"""Convenience routines for easily plotting with VisClaw.""" + +from __future__ import absolute_import +import os +import sys +import types +import six + +def plot(setplot="./sill.py", outdir="./_output", plotdir=None, htmlplot=True, + iplot=False, file_format='ascii', **plot_kargs): + r"""setplot can be a function or a path to a file.""" + + # Construct a plot directory if not provided + if plotdir is None: + try: + plotdir = os.path.join(outdir,"./_plots") + except AttributeError: + plotdir = os.path.join(os.getcwd(),"_plots") + + if htmlplot or iplot: + setplot_func = None + # No setplot specified, try to use a local file + if setplot is None: + # Grab and import the setplot function + local_setplot_path = os.path.join(os.getcwd(),'setplot.py') + if os.path.exists(local_setplot_path): + setplot = local_setplot_path + + # Fetch setplot function depending on type of setplot + if isinstance(setplot, types.FunctionType): + # setplot points to a function + setplot_func = lambda plotdata:setplot(plotdata, **plot_kargs) + + elif isinstance(setplot, types.ModuleType): + # setplot points to a module + setplot_func = lambda plotdata:setplot.setplot(plotdata, **plot_kargs) + + elif isinstance(setplot, six.string_types): + # setplot contains a path to a module + path = os.path.abspath(os.path.expandvars(os.path.expanduser(setplot))) + setplot_module_dir = os.path.dirname(path) + setplot_module_name = os.path.splitext(os.path.basename(setplot))[0] + sys.path.insert(0,setplot_module_dir) + setplot_module = __import__(setplot_module_name) + setplot_func = lambda plotdata:setplot_module.setplot(plotdata, **plot_kargs) + + if not isinstance(setplot_func, types.FunctionType): + # Everything else has failed, use default setplot + import clawpack.visclaw.setplot_default as setplot_module + setplot_func = setplot_module.setplot + + # Interactive plotting + if iplot: + from clawpack.visclaw import Iplotclaw + + ip = Iplotclaw.Iplotclaw(setplot=setplot_func, outdir=outdir) + ip.plotdata.format = file_format + + ip.plotloop() + + # Static HTML plotting + if htmlplot: + from clawpack.visclaw import plotclaw + plotclaw.plotclaw(outdir, plotdir, format=file_format, + setplot=setplot_func) + + +# These now just point to the above more generic function +def interactive_plot(outdir='./_output', file_format='ascii', setplot=None): + """Convenience function for launching an interactive plotting session.""" + plot(setplot, outdir=outdir, file_format=file_format, iplot=True, + htmlplot=False) + + +def html_plot(outdir='./_output', file_format='ascii', setplot=None): + """Convenience function for creating html page with plots.""" + plot(setplot, outdir=outdir, file_format=file_format, htmlplot=True, + iplot=False) + +plot() diff --git a/src/pyclaw/classic/solver.py b/src/pyclaw/classic/solver.py index 8b5634570..6f412065c 100644 --- a/src/pyclaw/classic/solver.py +++ b/src/pyclaw/classic/solver.py @@ -1,8 +1,8 @@ r""" Module containing the classic Clawpack solvers. -This module contains the pure and wrapped classic clawpack solvers. All -clawpack solvers inherit from the :class:`ClawSolver` superclass which in turn +This module contains the pure and wrapped classic clawpack solvers. All +clawpack solvers inherit from the :class:`ClawSolver` superclass which in turn inherits from the :class:`~pyclaw.solver.Solver` superclass. These are both pure virtual classes; the only solver classes that should be instantiated are the dimension-specific ones, :class:`ClawSolver1D` and :class:`ClawSolver2D`. @@ -21,38 +21,38 @@ class ClawSolver(Solver): r""" Generic classic Clawpack solver - + All Clawpack solvers inherit from this base class. - - .. attribute:: mthlim - + + .. attribute:: mthlim + Limiter(s) to be used. Specified either as one value or a list. If one value, the specified limiter is used for all wave families. If a list, the specified values indicate which limiter to apply to each wave family. Take a look at pyclaw.limiters.tvd for an enumeration. ``Default = limiters.tvd.minmod`` - + .. attribute:: order - + Order of the solver, either 1 for first order (i.e., Godunov's method) or 2 for second order (Lax-Wendroff-LeVeque). ``Default = 2`` - + .. attribute:: source_split - - Which source splitting method to use: 1 for first + + Which source splitting method to use: 1 for first order Godunov splitting and 2 for second order Strang splitting. ``Default = 1`` - + .. attribute:: fwave - - Whether to split the flux jump (rather than the jump in Q) into waves; - requires that the Riemann solver performs the splitting. + + Whether to split the flux jump (rather than the jump in Q) into waves; + requires that the Riemann solver performs the splitting. ``Default = False`` - + .. attribute:: step_source - - Handle for function that evaluates the source term. + + Handle for function that evaluates the source term. The required signature for this function is: def step_source(solver,state,dt) @@ -61,14 +61,14 @@ def step_source(solver,state,dt) Specifies whether to use wrapped Fortran routines ('Fortran') or pure Python ('Python'). ``Default = 'Fortran'``. - + .. attribute:: verbosity The level of detail of logged messages from the Fortran solver. ``Default = 0``. """ - + # ========== Generic Init Routine ======================================== def __init__(self,riemann_solver=None,claw_package=None): r""" @@ -93,7 +93,7 @@ def __init__(self,riemann_solver=None,claw_package=None): # Call general initialization function super(ClawSolver,self).__init__(riemann_solver,claw_package) - + # ========== Time stepping routines ====================================== def step(self,solution,take_one_step,tstart,tend): r""" @@ -102,15 +102,15 @@ def step(self,solution,take_one_step,tstart,tend): The elements of the algorithm for taking one step are: 1. Pick a step size as specified by the base solver attribute :func:`get_dt` - - 2. A half step on the source term :func:`step_source` if Strang splitting is + + 2. A half step on the source term :func:`step_source` if Strang splitting is being used (:attr:`source_split` = 2) - + 3. A step on the homogeneous problem :math:`q_t + f(q)_x = 0` is taken - + 4. A second half step or a full step is taken on the source term - :func:`step_source` depending on whether Strang splitting was used - (:attr:`source_split` = 2) or Godunov splitting + :func:`step_source` depending on whether Strang splitting was used + (:attr:`source_split` = 2) or Godunov splitting (:attr:`source_split` = 1) This routine is called from the method evolve_to_time defined in the @@ -118,8 +118,8 @@ def step(self,solution,take_one_step,tstart,tend): :Input: - *solution* - (:class:`~pyclaw.solution.Solution`) solution to be evolved - - :Output: + + :Output: - (bool) - True if full step succeeded, False otherwise """ self.get_dt(solution.t,tstart,tend,take_one_step) @@ -130,7 +130,7 @@ def step(self,solution,take_one_step,tstart,tend): self.step_hyperbolic(solution) - # Check here if the CFL condition is satisfied. + # Check here if the CFL condition is satisfied. # If not, return # immediately to evolve_to_time and let it deal with # picking a new step size (dt). if self.cfl.get_cached_max() >= self.cfl_max: @@ -144,7 +144,7 @@ def step(self,solution,take_one_step,tstart,tend): # Godunov Splitting if self.source_split == 1: self.step_source(self,solution.states[0],self.dt) - + return True def _check_cfl_settings(self): @@ -156,14 +156,14 @@ def _allocate_workspace(self,solution): def step_hyperbolic(self,solution): r""" Take one homogeneous step on the solution. - + This is a dummy routine and must be overridden. """ raise Exception("Dummy routine, please override!") def _set_mthlim(self): r""" - Convenience routine to convert users limiter specification to + Convenience routine to convert users limiter specification to the format understood by the Fortran code (i.e., a list of length num_waves). """ self._mthlim = self.limiters @@ -171,7 +171,7 @@ def _set_mthlim(self): if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves if len(self._mthlim)!=self.num_waves: raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves') - + def _set_method(self,state): r""" Set values of the solver._method array required by the Fortran code. @@ -251,13 +251,13 @@ def __del__(self): class ClawSolver1D(ClawSolver): r""" Clawpack evolution routine in 1D - - This class represents the 1d clawpack solver on a single grid. Note that - there are routines here for interfacing with the fortran time stepping - routines and the Python time stepping routines. The ones used are - dependent on the argument given to the initialization of the solver + + This class represents the 1d clawpack solver on a single grid. Note that + there are routines here for interfacing with the fortran time stepping + routines and the Python time stepping routines. The ones used are + dependent on the argument given to the initialization of the solver (defaults to python). - + """ __doc__ += add_parent_doc(ClawSolver) @@ -268,9 +268,9 @@ def __init__(self, riemann_solver=None, claw_package=None): Output: - (:class:`ClawSolver1D`) - Initialized 1d clawpack solver - + See :class:`ClawSolver1D` for more info. - """ + """ self.num_dim = 1 self.reflect_index = [1] @@ -283,7 +283,7 @@ def step_hyperbolic(self,solution): Take one time step on the homogeneous hyperbolic system. :Input: - - *solution* - (:class:`~pyclaw.solution.Solution`) Solution that + - *solution* - (:class:`~pyclaw.solution.Solution`) Solution that will be evolved """ import numpy as np @@ -292,32 +292,41 @@ def step_hyperbolic(self,solution): grid = state.grid self._apply_bcs(state) - + num_eqn,num_ghost = state.num_eqn,self.num_ghost - + if(self.kernel_language == 'Fortran'): mx = grid.num_cells[0] dx,dt = grid.delta[0],self.dt dtdx = np.zeros( (mx+2*num_ghost) ) + dt/dx rp1 = self.rp.rp1._cpointer - + self.qbc,cfl = self.fmod.step1(num_ghost,mx,self.qbc,self.auxbc,dx,dt,self._method,self._mthlim,self.fwave,rp1) - + elif(self.kernel_language == 'Python'): - + q = self.qbc aux = self.auxbc # Limiter to use in the pth family - limiter = np.array(self._mthlim,ndmin=1) - + limiter = np.array(self._mthlim,ndmin=1) + #print("num_ghost",self.num_ghost) dtdx = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) + #print(dtdx.shape) # Find local value for dt/dx if state.index_capa>=0: - dtdx = self.dt / (grid.delta[0] * state.aux[state.index_capa,:]) + # print(aux) + nw = 25 + cells_number = 100 + alpha = 0.25 + xpxc = (cells_number-1) * 1.0 / (cells_number) + aux[1,nw+2-1] = alpha *xpxc + aux[1,nw+2] = 1-alpha * xpxc + dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) #state.aux[state.index_capa,:]) + # print(dtdx.shape) else: dtdx += self.dt/grid.delta[0] - + # Solve Riemann problem at each interface q_l=q[:,:-1] q_r=q[:,1:] @@ -328,7 +337,7 @@ def step_hyperbolic(self,solution): aux_l = None aux_r = None wave,s,amdq,apdq = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - + # Update loop limits, these are the limits for the Riemann solver # locations, which then update a grid cell value # We include the Riemann problem just outside of the grid so we can @@ -337,14 +346,16 @@ def step_hyperbolic(self,solution): # | LL | | | | ... | | | UL | | # | | - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - + LL = self.num_ghost - 1 # 1 + # print("LL",LL) + UL = self.num_ghost + grid.num_cells[0] + 1 #35 + # print("UL",UL) + # print("dtdx",dtdx,"apdq",apdq.shape,"q",q.shape) # Update q for Godunov update for m in range(num_eqn): q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] - + # Compute maximum wave speed cfl = 0.0 for mw in range(wave.shape[1]): @@ -356,7 +367,7 @@ def step_hyperbolic(self,solution): if self.order == 2: # Initialize flux corrections f = np.zeros( (num_eqn,grid.num_cells[0] + 2*self.num_ghost) ) - + # Apply Limiters to waves if (limiter > 0).any(): wave = tvd.limit(state.num_eqn,wave,s,limiter,dtdx) @@ -379,7 +390,7 @@ def step_hyperbolic(self,solution): # Update q by differencing correction fluxes for m in range(num_eqn): - q[m,LL:UL-1] -= dtdx[LL:UL-1] * (f[m,LL+1:UL] - f[m,LL:UL-1]) + q[m,LL:UL-1] -= dtdx[LL:UL-1] * (f[m,LL+1:UL] - f[m,LL:UL-1]) else: raise Exception("Unrecognized kernel_language; choose 'Fortran' or 'Python'") @@ -387,7 +398,7 @@ def step_hyperbolic(self,solution): state.set_q_from_qbc(num_ghost,self.qbc) if state.num_aux > 0: state.set_aux_from_auxbc(num_ghost,self.auxbc) - + # ============================================================================ # ClawPack 2d Solver Class @@ -401,9 +412,9 @@ class ClawSolver2D(ClawSolver): In addition to the attributes of ClawSolver1D, ClawSolver2D also has the following options: - + .. attribute:: dimensional_split - + If True, use dimensional splitting (Godunov splitting). Dimensional splitting with Strang splitting is not supported at present but could easily be enabled if necessary. @@ -411,14 +422,14 @@ class ClawSolver2D(ClawSolver): transverse Riemann solves. .. attribute:: transverse_waves - + If dimensional_split is True, this option has no effect. If dimensional_split is False, then transverse_waves should be one of the following values: ClawSolver2D.no_trans: Transverse Riemann solver not used. The stable CFL for this algorithm is 0.5. Not recommended. - + ClawSolver2D.trans_inc: Transverse increment waves are computed and propagated. @@ -429,7 +440,7 @@ class ClawSolver2D(ClawSolver): """ __doc__ += add_parent_doc(ClawSolver) - + no_trans = 0 trans_inc = 1 trans_cor = 2 @@ -437,9 +448,9 @@ class ClawSolver2D(ClawSolver): def __init__(self,riemann_solver=None, claw_package=None): r""" Create 2d Clawpack solver - + See :class:`ClawSolver2D` for more info. - """ + """ self.dimensional_split = True self.transverse_waves = self.trans_inc @@ -511,10 +522,10 @@ def step_hyperbolic(self,solution): dx,dy = grid.delta mx,my = grid.num_cells maxm = max(mx,my) - + self._apply_bcs(state) qold = self.qbc.copy('F') - + rpn2 = self.rp.rpn2._cpointer if (self.dimensional_split) or (self.transverse_waves==0): @@ -562,9 +573,9 @@ class ClawSolver3D(ClawSolver): In addition to the attributes of ClawSolver, ClawSolver3D also has the following options: - + .. attribute:: dimensional_split - + If True, use dimensional splitting (Godunov splitting). Dimensional splitting with Strang splitting is not supported at present but could easily be enabled if necessary. @@ -572,14 +583,14 @@ class ClawSolver3D(ClawSolver): transverse Riemann solves. .. attribute:: transverse_waves - + If dimensional_split is True, this option has no effect. If dim_plit is False, then transverse_waves should be one of the following values: ClawSolver3D.no_trans: Transverse Riemann solver not used. The stable CFL for this algorithm is 0.5. Not recommended. - + ClawSolver3D.trans_inc: Transverse increment waves are computed and propagated. @@ -599,9 +610,9 @@ class ClawSolver3D(ClawSolver): def __init__(self, riemann_solver=None, claw_package=None): r""" Create 3d Clawpack solver - + See :class:`ClawSolver3D` for more info. - """ + """ # Add the functions as required attributes self.dimensional_split = True self.transverse_waves = self.trans_cor @@ -616,7 +627,7 @@ def __init__(self, riemann_solver=None, claw_package=None): super(ClawSolver3D,self).__init__(riemann_solver, claw_package) - # ========== Setup routine ============================= + # ========== Setup routine ============================= def _allocate_workspace(self,solution): r""" Allocate auxN and work arrays for use in Fortran subroutines. @@ -661,11 +672,11 @@ def step_hyperbolic(self,solution): dx,dy,dz = grid.delta mx,my,mz = grid.num_cells maxm = max(mx,my,mz) - + self._apply_bcs(state) qnew = self.qbc qold = qnew.copy('F') - + rpn3 = self.rp.rpn3._cpointer if (self.dimensional_split) or (self.transverse_waves==0): From 174ff78ba97ffd7101b8d0b3119db07752c848ec Mon Sep 17 00:00:00 2001 From: cr2940 Date: Tue, 8 Oct 2019 12:08:12 -0400 Subject: [PATCH 02/47] added nonuniform example --- examples/advection_1d/advection_1d.py | 48 ++++--- examples/advection_1d/advection_1d_nonunif.py | 125 ++++++++++++++++++ 2 files changed, 156 insertions(+), 17 deletions(-) create mode 100644 examples/advection_1d/advection_1d_nonunif.py diff --git a/examples/advection_1d/advection_1d.py b/examples/advection_1d/advection_1d.py index 533505165..10d661c2f 100755 --- a/examples/advection_1d/advection_1d.py +++ b/examples/advection_1d/advection_1d.py @@ -7,7 +7,7 @@ Solve the linear advection equation: -.. math:: +.. math:: q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. @@ -19,9 +19,13 @@ from __future__ import absolute_import import numpy as np from clawpack import riemann +from clawpack.pyclaw.plot import plot -def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, - time_integrator='SSP104', outdir='./_output'): + + + +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, + time_integrator='SSP104', outdir='./_output_hb'): if use_petsc: import clawpack.petclaw as pyclaw @@ -32,7 +36,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi riemann_solver = riemann.advection_1D elif kernel_language == 'Python': riemann_solver = riemann.advection_1D_py.advection_1D - + if solver_type=='classic': solver = pyclaw.ClawSolver1D(riemann_solver) elif solver_type=='sharpclaw': @@ -45,13 +49,18 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language - + solver.order=1 + solver.limiters = None + solver.num_eqn=1 + solver.num_waves=1 solver.bc_lower[0] = pyclaw.BC.periodic solver.bc_upper[0] = pyclaw.BC.periodic + solver.aux_bc_lower[0] = pyclaw.BC.periodic + solver.aux_bc_upper[0] = pyclaw.BC.periodic x = pyclaw.Dimension(0.0,1.0,nx,name='x') domain = pyclaw.Domain(x) - state = pyclaw.State(domain,solver.num_eqn) + state = pyclaw.State(domain,1) state.problem_data['u'] = 1. # Advection velocity @@ -65,20 +74,23 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi claw.solution = pyclaw.Solution(state,domain) claw.solver = solver - if outdir is not None: - claw.outdir = outdir - else: + claw.tfinal =1.0 + claw.outdir = outdir + if outdir is None: claw.output_format = None - claw.tfinal =1.0 + claw.run() + + claw.setplot = setplot - return claw + + # return claw def setplot(plotdata): - """ + """ Plot solution using VisClaw. - """ + """ plotdata.clearfigures() # clear any old figures,axes,items data plotfigure = plotdata.new_plotfigure(name='q', figno=1) @@ -94,10 +106,12 @@ def setplot(plotdata): plotitem.plotstyle = '-o' plotitem.color = 'b' plotitem.kwargs = {'linewidth':2,'markersize':5} - + return plotdata - +#plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) + if __name__=="__main__": - from clawpack.pyclaw.util import run_app_from_main - output = run_app_from_main(setup,setplot) + #from clawpack.pyclaw.util import run_app_from_main + #output = run_app_from_main(setup,setplot) + setup() diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py new file mode 100644 index 000000000..9d7c83724 --- /dev/null +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python +# encoding: utf-8 + +r""" +One-dimensional advection +========================= + +Solve the linear advection equation: + +.. math:: + q_t + u q_x = 0. + +Here q is the density of some conserved quantity and u is the velocity. + +The initial condition is a Gaussian and the boundary conditions are periodic. +The final solution is identical to the initial data because the wave has +crossed the domain exactly once. +""" +from __future__ import absolute_import +import numpy as np +from clawpack import riemann +from clawpack.pyclaw.plot import plot + + + + +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, + time_integrator='SSP104', outdir='./_output_hb'): + + if use_petsc: + import clawpack.petclaw as pyclaw + else: + from clawpack import pyclaw + + if kernel_language == 'Fortran': + riemann_solver = riemann.advection_1D + elif kernel_language == 'Python': + riemann_solver = riemann.advection_1D_py.advection_1D + + if solver_type=='classic': + solver = pyclaw.ClawSolver1D(riemann_solver) + elif solver_type=='sharpclaw': + solver = pyclaw.SharpClawSolver1D(riemann_solver) + solver.weno_order = weno_order + solver.time_integrator = time_integrator + if time_integrator == 'SSPLMMk3': + solver.lmm_steps = 5 + solver.check_lmm_cond = True + else: raise Exception('Unrecognized value of solver_type.') + + solver.kernel_language = kernel_language + solver.order=1 + solver.limiters = None + solver.num_eqn=1 + solver.num_waves=1 + solver.bc_lower[0] = pyclaw.BC.periodic + solver.bc_upper[0] = pyclaw.BC.periodic + solver.aux_bc_lower[0] = pyclaw.BC.periodic + solver.aux_bc_upper[0] = pyclaw.BC.periodic + + x = pyclaw.Dimension(0.0,1.0,nx,name='x') + domain = pyclaw.Domain(x) + state = pyclaw.State(domain,1) + + state.problem_data['u'] = 1. # Advection velocity + nw = 50 # location of nouniformity + alpha = 0.1 # ratio of nonuniform small cell to standard cell + + state.index_capa = 0 + xpxc = nx * 1.0 / (nx-1) + state.aux = np.zeros((1,nx)) + state.aux[0, :] = xpxc + state.aux[0, nw-1] = alpha * xpxc + state.aux[0, nw] = (1-alpha) * xpxc + + # Initial data + xc = state.grid.x.centers + beta = 100; gamma = 0; x0 = 0.75 + state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) + + claw = pyclaw.Controller() + claw.keep_copy = True + claw.solution = pyclaw.Solution(state,domain) + claw.solver = solver + + claw.tfinal =1.0 + claw.outdir = outdir + if outdir is None: + claw.output_format = None + + claw.run() + + + claw.setplot = setplot + +# plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) + + # return claw + +def setplot(plotdata): + """ + Plot solution using VisClaw. + """ + plotdata.clearfigures() # clear any old figures,axes,items data + + plotfigure = plotdata.new_plotfigure(name='q', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.ylimits = [-.2,1.0] + plotaxes.title = 'q' + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = '-o' + plotitem.color = 'b' + plotitem.kwargs = {'linewidth':2,'markersize':5} + + return plotdata + +if __name__=="__main__": + #from clawpack.pyclaw.util import run_app_from_main + #output = run_app_from_main(setup,setplot) + setup() From 465872ca00d06cbb627bfbce4d0f410714357bcc Mon Sep 17 00:00:00 2001 From: cr2940 Date: Tue, 8 Oct 2019 12:08:59 -0400 Subject: [PATCH 03/47] added new nonuniform example --- examples/shallow_1d/dam_break.py | 4 + src/pyclaw/classic/solver.py | 26 ++---- src/pyclaw/limiters/tvd.py | 154 +++++++++++++++---------------- 3 files changed, 87 insertions(+), 97 deletions(-) diff --git a/examples/shallow_1d/dam_break.py b/examples/shallow_1d/dam_break.py index 835040b0e..cb09ad54a 100755 --- a/examples/shallow_1d/dam_break.py +++ b/examples/shallow_1d/dam_break.py @@ -74,6 +74,7 @@ def setup(use_petsc=False,kernel_language='Fortran',outdir='./_output',solver_ty ur = 0. state.q[depth,:] = hl * (xc <= x0) + hr * (xc > x0) state.q[momentum,:] = hl*ul * (xc <= x0) + hr*ur * (xc > x0) + elif IC=='2-shock': hl = 1. ul = 1. @@ -81,10 +82,13 @@ def setup(use_petsc=False,kernel_language='Fortran',outdir='./_output',solver_ty ur = -1. state.q[depth,:] = hl * (xc <= x0) + hr * (xc > x0) state.q[momentum,:] = hl*ul * (xc <= x0) + hr*ur * (xc > x0) + elif IC=='perturbation': eps=0.1 state.q[depth,:] = 1.0 + eps*np.exp(-(xc-x0)**2/0.5) state.q[momentum,:] = 0. + state.aux[0,:] = np.zeros((1,mx)) + state.aux[0,int(mx/2):] +=0.5 claw = pyclaw.Controller() claw.keep_copy = True diff --git a/src/pyclaw/classic/solver.py b/src/pyclaw/classic/solver.py index 6f412065c..cbceb2827 100644 --- a/src/pyclaw/classic/solver.py +++ b/src/pyclaw/classic/solver.py @@ -79,7 +79,7 @@ def __init__(self,riemann_solver=None,claw_package=None): """ self.num_ghost = 2 self.limiters = tvd.minmod - self.order = 2 + self.order = 1 self.source_split = 1 self.fwave = False self.step_source = None @@ -276,7 +276,6 @@ def __init__(self, riemann_solver=None, claw_package=None): super(ClawSolver1D,self).__init__(riemann_solver, claw_package) - # ========== Homogeneous Step ===================================== def step_hyperbolic(self,solution): r""" @@ -290,6 +289,7 @@ def step_hyperbolic(self,solution): state = solution.states[0] grid = state.grid + #print(grid) self._apply_bcs(state) @@ -307,23 +307,12 @@ def step_hyperbolic(self,solution): q = self.qbc aux = self.auxbc - # Limiter to use in the pth family limiter = np.array(self._mthlim,ndmin=1) - #print("num_ghost",self.num_ghost) dtdx = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - #print(dtdx.shape) # Find local value for dt/dx if state.index_capa>=0: - # print(aux) - nw = 25 - cells_number = 100 - alpha = 0.25 - xpxc = (cells_number-1) * 1.0 / (cells_number) - aux[1,nw+2-1] = alpha *xpxc - aux[1,nw+2] = 1-alpha * xpxc dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) #state.aux[state.index_capa,:]) - # print(dtdx.shape) else: dtdx += self.dt/grid.delta[0] @@ -346,16 +335,14 @@ def step_hyperbolic(self,solution): # | LL | | | | ... | | | UL | | # | | - LL = self.num_ghost - 1 # 1 - # print("LL",LL) - UL = self.num_ghost + grid.num_cells[0] + 1 #35 - # print("UL",UL) - # print("dtdx",dtdx,"apdq",apdq.shape,"q",q.shape) - # Update q for Godunov update + LL = self.num_ghost - 1 + UL = self.num_ghost + grid.num_cells[0] + 1 + for m in range(num_eqn): q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] + # Compute maximum wave speed cfl = 0.0 for mw in range(wave.shape[1]): @@ -399,7 +386,6 @@ def step_hyperbolic(self,solution): if state.num_aux > 0: state.set_aux_from_auxbc(num_ghost,self.auxbc) - # ============================================================================ # ClawPack 2d Solver Class # ============================================================================ diff --git a/src/pyclaw/limiters/tvd.py b/src/pyclaw/limiters/tvd.py index 5195cacb8..af49ba11f 100644 --- a/src/pyclaw/limiters/tvd.py +++ b/src/pyclaw/limiters/tvd.py @@ -22,51 +22,51 @@ 6. Frommm - :math:`1/2 (1 + r)` 7. Albada 2 - :math:`(r^2 + r) / (1 + r^2)` 8. Albada 3 - :math:`1/2 (1+r) (1 - (|1-r|^3) / (1+|r|^3))` -9. van Leer with Klein sharpening, k=2 - +9. van Leer with Klein sharpening, k=2 - :func:`van_leer_klein_sharpening_limiter` CFL Dependent Limiters '''''''''''''''''''''' 10. Roe's linear third order scheme - :math:`1 + (r-1) (1 + cfl) / 3` -11. Arora-Roe (= limited version of the linear third order scheme) - +11. Arora-Roe (= limited version of the linear third order scheme) - :func:`arora_roe` -12. Theta Limiter, theta=0.95 (safety on nonlinear waves) - +12. Theta Limiter, theta=0.95 (safety on nonlinear waves) - :func:`theta_limiter` 13. Theta Limiter, theta=0.75 - :func:`theta_limiter` 14. Theta Limiter, theta=0.5 - :func:`theta_limiter` 15. CFL-Superbee (Roe's Ultrabee) - :func:`cfl_superbee` -16. CFL-Superbee (Roe's Ultrabee) with theta=0.95 (nonlinear waves) - +16. CFL-Superbee (Roe's Ultrabee) with theta=0.95 (nonlinear waves) - :func:`cfl_superbee_theta` 17. beta=2/3 limiter - :func:`beta_limiter` 18. beta=2/3 limiter with theta=0.95 (nonlinear waves) - :func:`beta_limiter` 19. Hyperbee - :func:`hyperbee_limiter` 20. SuperPower - :func:`superpower_limiter` 21. Cada-Torrilhon modified - :func:`cada_torrilhon_limiter` -22. Cada-Torrilhon modified, version for nonlinear waves - +22. Cada-Torrilhon modified, version for nonlinear waves - :func:`cada_torrilhon_limiter_nonlinear` 23. upper bound limiter (1st order) - :func:`upper_bound_limiter` - + All limiters have the same function call signature: :Input: - - *r* - (ndarray(:)) + - *r* - (ndarray(:)) - *cfl* - (ndarray(:)) Local CFL number - + :Output: - - (ndarray(:)) - + - (ndarray(:)) - -Newer limiters are based on work done by Friedemann Kemm [kemm_2009]_, paper +Newer limiters are based on work done by Friedemann Kemm [kemm_2009]_, paper in review. - + :Authors: Kyle Mandli and Randy LeVeque (2008-08-21) Initial version - + Kyle Mandli (2009-07-05) Added CFL depdendent limiters """ # ============================================================================ # Copyright (C) 2008 Kyle T. Mandli # Copyrigth (C) 2009 Randall J. LeVeque # -# Distributed under the terms of the Berkeley Software Distribution (BSD) +# Distributed under the terms of the Berkeley Software Distribution (BSD) # license # http://www.opensource.org/licenses/ # ============================================================================ @@ -85,25 +85,25 @@ def limit(num_eqn,wave,s,limiter,dtdx): Apply a limiter to the waves Function that limits the given waves using the methods contained - in limiter. This is the vectorized version of the function acting on a + in limiter. This is the vectorized version of the function acting on a row of waves at a time. - + :Input: - *wave* - (ndarray(:,num_eqn,num_waves)) The waves at each interface - *s* - (ndarray(:,num_waves)) Speeds for each wave - - *limiter* - (``int`` list) Array of type ``int`` determining which + - *limiter* - (``int`` list) Array of type ``int`` determining which limiter to use - - *dtdx* - (ndarray(:)) :math:`\Delta t / \Delta x` ratio, used for CFL + - *dtdx* - (ndarray(:)) :math:`\Delta t / \Delta x` ratio, used for CFL dependent limiters - + :Output: - (ndarray(:,num_eqn,num_waves)) - Returns the limited waves :Version: 1.1 (2009-07-05) """ - + # wave_norm2 is the sum of the squares along the num_eqn axis, - # so the norm of the cell i for wave number j is addressed + # so the norm of the cell i for wave number j is addressed # as wave_norm2[i,j] wave_norm2 = np.sum(np.square(wave),axis=0) wave_zero_mask = np.array((wave_norm2 == 0), dtype=float) @@ -126,13 +126,13 @@ def limit(num_eqn,wave,s,limiter,dtdx): # Fill the rest of the array r.fill_value = 0 r = r.filled() - + for mw in range(wave.shape[1]): # skip waves that are marked as not needing a limiter limit_func = limiter_functions.get(limiter[mw]) if limit_func is not None: for m in range(num_eqn): - cfl = np.abs(s[mw,1:-1]*(dtdx[1:-2]*spos[mw,:] + cfl = np.abs(s[mw,1:-1]*(dtdx[1:-2]*spos[mw,:] + (1-spos[mw,:])*dtdx[2:-1])) wlimitr = limit_func(r[mw,:],cfl) wave[m,mw,1:-1] = wave[m,mw,1:-1]*wave_zero_mask[mw,1:-1] \ @@ -146,13 +146,13 @@ def minmod_limiter(r,cfl): """ a = np.ones((2,len(r))) b = np.zeros((2,len(r))) - + a[1,:] = r - + b[1,:] = np.minimum(a[0],a[1]) - + return np.maximum(b[0],b[1]) - + def superbee_limiter(r,cfl): r""" @@ -163,14 +163,14 @@ def superbee_limiter(r,cfl): c = np.zeros((3,len(r))) a[1,:] = 2.0*r - + b[1,:] = r c[1,:] = np.minimum(a[0], a[1]) c[2,:] = np.minimum(b[0], b[1]) return np.max(c,axis=0) - + def mc_limiter(r,cfl): r""" MC vectorized limiter @@ -181,23 +181,23 @@ def mc_limiter(r,cfl): a[0,:] = (1.0 + r) / 2.0 a[1,:] = 2 a[2,:] = 2.0 * r - + b[1,:] = np.min(a,axis=0) return np.maximum(b[0], b[1]) - + def van_leer_klein_sharpening_limiter(r,cfl): r""" van Leer with Klein sharpening, k=2 """ a = np.ones((2,len(r))) * 1.e-5 a[0,:] = r - + rcorr = np.maximum(a[0],a[1]) a[1,:] = 1/rcorr sharg = np.minimum(a[0],a[1]) sharp = 1.0 + sharg * (1.0 - sharg) * (1.0 - sharg**2) - + return (r+np.abs(r)) / (1.0 + np.abs(r)) * sharp def arora_roe(r,cfl): @@ -205,31 +205,31 @@ def arora_roe(r,cfl): Arora-Roe limiter, limited version of the linear third order scheme """ caut = 0.99 - + a = np.empty((3,len(r))) b = np.zeros((2,len(r))) - + s1 = (caut * 2.0 / cfl) s2 = (1.0 + cfl) / 3.0 phimax = caut * 2.0 / (1.0 - cfl) - + a[0,:] = s1 * r a[1,:] = 1.0 + s2 * (r - 1.0) a[2,:] = phimax b[1,:] = np.min(a,axis=0) - + return np.maximum(b[0],b[1]) def theta_limiter(r,cfl,theta=0.95): r""" Theta limiter - + Additional Input: - *theta* = """ a = np.empty((2,len(r))) b = np.empty((3,len(r))) - + a[0,:] = 0.001 a[1,:] = cfl cfmod1 = np.max(a,axis=0) @@ -238,50 +238,50 @@ def theta_limiter(r,cfl,theta=0.95): s1 = 2.0 / cfmod1 s2 = (1.0 + cfl) / 3.0 phimax = 2.0 / (1.0 - cfmod2) - + a[0,:] = (1.0 - theta) * s1 a[1,:] = 1.0 + s2 * (r - 1.0) left = np.maximum(a[0],a[1]) a[0,:] = (1.0 - theta) * phimax * r a[1,:] = theta * s1 * r middle = np.maximum(a[0],a[1]) - + b[0,:] = left b[1,:] = middle b[2,:] = theta*phimax - + return np.min(b,axis=0) - + def cfl_superbee(r,cfl): r""" CFL-Superbee (Roe's Ultrabee) without theta parameter """ a = np.empty((2,len(r))) b = np.zeros((3,len(r))) - + a[0,:] = 0.001 a[1,:] = cfl cfmod1 = np.maximum(a[0],a[1]) a[0,:] = 0.999 cfmod2 = np.minimum(a[0],a[1]) - + a[0,:] = 1.0 a[1,:] = 2.0 * r / cfmod1 b[1,:] = np.minimum(a[0],a[1]) a[0,:] = 2.0/(1-cfmod2) a[1,:] = r b[2,:] = np.minimum(a[0],a[1]) - + return np.max(b,axis=0) - + def cfl_superbee_theta(r,cfl,theta=0.95): r""" CFL-Superbee (Roe's Ultrabee) with theta parameter """ a = np.empty((2,len(r))) b = np.zeros((2,len(r))) - + a[0,:] = 0.001 a[1,:] = cfl cfmod1 = np.maximum(a[0],a[1]) @@ -295,7 +295,7 @@ def cfl_superbee_theta(r,cfl,theta=0.95): a[1,:] = phimax b[1,:] = np.minimum(a[0],a[1]) ultra = np.maximum(b[0],b[1]) - + a[0,:] = ultra b[0,:] = 1.0 b[1,:] = r @@ -305,49 +305,49 @@ def cfl_superbee_theta(r,cfl,theta=0.95): def beta_limiter(r,cfl,theta=0.95,beta=0.66666666666666666): r""" Modification of CFL Superbee limiter with theta and beta parameters - + Additional Input: - *theta* - *beta* """ a = np.empty((2,len(r))) b = np.zeros((2,len(r))) - + a[0,:] = 0.001 a[1,:] = cfl cfmod1 = np.maximum(a[0],a[1]) a[0,:] = 0.999 cfmod2 = np.minimum(a[0],a[1]) - + s1 = theta * 2.0 / cfmod1 s2 = (1.0 + cfl) / 3.0 phimax = theta * 2.0 / (1.0 - cfmod2) - + a[0,:] = s1*r a[1,:] = phimax b[1,:] = np.min(a) ultra = np.max(b) - + a[0,:] = 1.0 + (s2 - beta/2.0) * (r-1.0) a[1,:] = 1.0 + (s2 + beta/2.0) * (r-1.0) b[0,:] = ultra b[1,:] = np.max(a) a[0,:] = 0.0 a[1,:] = np.min(b) - + return np.max(a) def hyperbee_limiter(r,cfl): r"""Hyperbee""" a = np.empty((2,len(r))) - + a[0,:] = 0.001 a[1,:] = cfl cfmod1 = np.maximum(a[0],a[1]) a[0,:] = 0.999 cfmod2 = np.minimum(a[0],a[1]) - + index1 = r < 0.0 index2 = np.abs(r-1.0) < 1.0e-6 index3 = np.abs(index1 + index2 - 1) @@ -356,34 +356,34 @@ def hyperbee_limiter(r,cfl): rmin = r-1.0 rdur = r/rmin return np.choose(master_index,[0.0,1.0, - 2.0 * rdur * (cfl * rmin + 1.0 - r**cfl) + 2.0 * rdur * (cfl * rmin + 1.0 - r**cfl) / (cfmod1 * (1.0 - cfmod2) * rmin)]) - + def superpower_limiter(r,cfl,caut=1.0): r""" SuperPower limiter - + Additional input: - *caut* = Limiter parameter """ s2 = (1.0 + cfl) / 3.0 s3 = 1.0 - s2 - - pp = (((r<=1.0) * np.abs(2.0/cfl) * caut) * 2.0 * s3 + + pp = (((r<=1.0) * np.abs(2.0/cfl) * caut) * 2.0 * s3 + (r > 1.0) * (np.abs(2.0)/(1.0 - cfl) * caut) * 2.0 * s2) - + rabs = np.abs(r) rfrac = np.abs((1.0 - rabs) / (1.0 + rabs)) signum = np.floor(0.5 * (1.0 + np.sign(r))) - + return signum * (s3 + s2 * r) * (1.0 - rfrac**pp) - + def cada_torrilhon_limiter(r,cfl,epsilon=1.0e-3): r""" Cada-Torrilhon modified - + Additional Input: - - *epsilon* = + - *epsilon* = """ a = np.ones((2,len(r))) * 0.95 b = np.empty((3,len(r))) @@ -392,7 +392,7 @@ def cada_torrilhon_limiter(r,cfl,epsilon=1.0e-3): cfl = np.min(a) a[1,:] = 0.05 cfl = np.max(a) - + # Multiply all parts except b[0,:] by (1.0 - epsilon) as well b[0,:] = 1.0 + (1+cfl) / 3.0 * (r - 1) b[1,:] = 2.0 * np.abs(r) / (cfl + epsilon) @@ -401,48 +401,48 @@ def cada_torrilhon_limiter(r,cfl,epsilon=1.0e-3): a[0,:] = np.min(b) a[1,:] = (-2.0 * (cfl**2 - 3.0 * cfl + 8.0) * (1.0-epsilon) / (np.abs(r) * (cfl**3 - cfl**2 - cfl + 1.0 + epsilon))) - + return np.max(a) - + def cada_torrilhon_limiter_nonlinear(r,cfl): r""" Cada-Torrilhon modified, version for nonlinear waves """ a = np.empty((3,len(r))) b = np.empty((2,len(r))) - + s2 = (1.0 + cfl) / 3.0 a[0,:] = 1.0 + s2 * (r - 1.0) a[1,:] = 2.0 * np.abs(r) / (0.6 + np.abs(r)) a[2,:] = 5.0 / np.abs(r) b[0,:] = np.min(a) b[1,:] = -3.0 / np.abs(r) - + return np.max(b) - + def upper_bound_limiter(r,cfl,theta=1.0): r""" Upper bound limiter (1st order) - + Additional Input: - *theta* = """ a = np.empty((2,len(r))) b = np.zeros((2,len(r))) - + a[0,:] = 0.001 a[1,:] = cfl cfmod1 = np.maximum(a[0],a[1]) a[0,:] = 0.999 cfmod2 = np.minimum(a[0],a[1]) - + s1 = theta * 2.0 / cfmod1 phimax = theta * 2.0 / (1.0 - cfmod2) - + a[0,:] = s1*r a[1,:] = phimax b[1,:] = np.min(a) - + return np.max(b) From dedd06b9de499ad09dfb3795793123a91540f661 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Tue, 8 Oct 2019 12:14:26 -0400 Subject: [PATCH 04/47] added nonunif example test file --- .../advection_1d/test_advection_nonunif.py | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 examples/advection_1d/test_advection_nonunif.py diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py new file mode 100644 index 000000000..f0e439a74 --- /dev/null +++ b/examples/advection_1d/test_advection_nonunif.py @@ -0,0 +1,52 @@ +from __future__ import absolute_import +def test_1d_advection(): + """test_1d_advection + + tests against expected classic, sharpclaw, and high-order weno results """ + + from . import advection_1d + + def verify_expected(expected): + """ given an expected value, returns a verification function """ + def advection_verify(claw): + from clawpack.pyclaw.util import check_diff + import numpy as np + + q0=claw.frames[0].state.get_q_global() + qfinal=claw.frames[claw.num_output_times].state.get_q_global() + + if q0 is not None and qfinal is not None: + dx=claw.solution.domain.grid.delta[0] + test = dx*np.linalg.norm(qfinal-q0,1) + return check_diff(expected, test, reltol=1e-4) + else: + return + return advection_verify + + from clawpack.pyclaw.util import gen_variants + + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(3.203924e-04), + kernel_languages=('Python','Fortran'), + solver_type='classic', outdir=None) + + sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected(1.163605e-05), + kernel_languages=('Python','Fortran'), + solver_type='sharpclaw',time_integrator='SSP104', outdir=None) + + sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected(1.500727e-05), + kernel_languages=('Python','Fortran'), + solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) + + weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(7.489618e-06), + kernel_languages=('Fortran',), solver_type='sharpclaw', + time_integrator='SSP104', weno_order=17, + outdir=None) + + from itertools import chain + for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests): + yield test + + +if __name__=="__main__": + import nose + nose.main() From 0bdee1268d530a15708f2281b6548f3df1ec5b9b Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 09:19:34 -0700 Subject: [PATCH 05/47] nonuniform advection example files --- examples/advection_1d/advection_1d_nonunif.py | 125 ++++++++++++++++++ .../advection_1d/test_advection_nonunif.py | 52 ++++++++ 2 files changed, 177 insertions(+) create mode 100644 examples/advection_1d/advection_1d_nonunif.py create mode 100644 examples/advection_1d/test_advection_nonunif.py diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py new file mode 100644 index 000000000..9d7c83724 --- /dev/null +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python +# encoding: utf-8 + +r""" +One-dimensional advection +========================= + +Solve the linear advection equation: + +.. math:: + q_t + u q_x = 0. + +Here q is the density of some conserved quantity and u is the velocity. + +The initial condition is a Gaussian and the boundary conditions are periodic. +The final solution is identical to the initial data because the wave has +crossed the domain exactly once. +""" +from __future__ import absolute_import +import numpy as np +from clawpack import riemann +from clawpack.pyclaw.plot import plot + + + + +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, + time_integrator='SSP104', outdir='./_output_hb'): + + if use_petsc: + import clawpack.petclaw as pyclaw + else: + from clawpack import pyclaw + + if kernel_language == 'Fortran': + riemann_solver = riemann.advection_1D + elif kernel_language == 'Python': + riemann_solver = riemann.advection_1D_py.advection_1D + + if solver_type=='classic': + solver = pyclaw.ClawSolver1D(riemann_solver) + elif solver_type=='sharpclaw': + solver = pyclaw.SharpClawSolver1D(riemann_solver) + solver.weno_order = weno_order + solver.time_integrator = time_integrator + if time_integrator == 'SSPLMMk3': + solver.lmm_steps = 5 + solver.check_lmm_cond = True + else: raise Exception('Unrecognized value of solver_type.') + + solver.kernel_language = kernel_language + solver.order=1 + solver.limiters = None + solver.num_eqn=1 + solver.num_waves=1 + solver.bc_lower[0] = pyclaw.BC.periodic + solver.bc_upper[0] = pyclaw.BC.periodic + solver.aux_bc_lower[0] = pyclaw.BC.periodic + solver.aux_bc_upper[0] = pyclaw.BC.periodic + + x = pyclaw.Dimension(0.0,1.0,nx,name='x') + domain = pyclaw.Domain(x) + state = pyclaw.State(domain,1) + + state.problem_data['u'] = 1. # Advection velocity + nw = 50 # location of nouniformity + alpha = 0.1 # ratio of nonuniform small cell to standard cell + + state.index_capa = 0 + xpxc = nx * 1.0 / (nx-1) + state.aux = np.zeros((1,nx)) + state.aux[0, :] = xpxc + state.aux[0, nw-1] = alpha * xpxc + state.aux[0, nw] = (1-alpha) * xpxc + + # Initial data + xc = state.grid.x.centers + beta = 100; gamma = 0; x0 = 0.75 + state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) + + claw = pyclaw.Controller() + claw.keep_copy = True + claw.solution = pyclaw.Solution(state,domain) + claw.solver = solver + + claw.tfinal =1.0 + claw.outdir = outdir + if outdir is None: + claw.output_format = None + + claw.run() + + + claw.setplot = setplot + +# plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) + + # return claw + +def setplot(plotdata): + """ + Plot solution using VisClaw. + """ + plotdata.clearfigures() # clear any old figures,axes,items data + + plotfigure = plotdata.new_plotfigure(name='q', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.ylimits = [-.2,1.0] + plotaxes.title = 'q' + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = '-o' + plotitem.color = 'b' + plotitem.kwargs = {'linewidth':2,'markersize':5} + + return plotdata + +if __name__=="__main__": + #from clawpack.pyclaw.util import run_app_from_main + #output = run_app_from_main(setup,setplot) + setup() diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py new file mode 100644 index 000000000..f0e439a74 --- /dev/null +++ b/examples/advection_1d/test_advection_nonunif.py @@ -0,0 +1,52 @@ +from __future__ import absolute_import +def test_1d_advection(): + """test_1d_advection + + tests against expected classic, sharpclaw, and high-order weno results """ + + from . import advection_1d + + def verify_expected(expected): + """ given an expected value, returns a verification function """ + def advection_verify(claw): + from clawpack.pyclaw.util import check_diff + import numpy as np + + q0=claw.frames[0].state.get_q_global() + qfinal=claw.frames[claw.num_output_times].state.get_q_global() + + if q0 is not None and qfinal is not None: + dx=claw.solution.domain.grid.delta[0] + test = dx*np.linalg.norm(qfinal-q0,1) + return check_diff(expected, test, reltol=1e-4) + else: + return + return advection_verify + + from clawpack.pyclaw.util import gen_variants + + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(3.203924e-04), + kernel_languages=('Python','Fortran'), + solver_type='classic', outdir=None) + + sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected(1.163605e-05), + kernel_languages=('Python','Fortran'), + solver_type='sharpclaw',time_integrator='SSP104', outdir=None) + + sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected(1.500727e-05), + kernel_languages=('Python','Fortran'), + solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) + + weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(7.489618e-06), + kernel_languages=('Fortran',), solver_type='sharpclaw', + time_integrator='SSP104', weno_order=17, + outdir=None) + + from itertools import chain + for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests): + yield test + + +if __name__=="__main__": + import nose + nose.main() From d3374ba1092c28caa1dc49de515826730e98007c Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 12:20:26 -0400 Subject: [PATCH 06/47] fix aux bug --- src/pyclaw/classic/solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyclaw/classic/solver.py b/src/pyclaw/classic/solver.py index 8b5634570..3f104e4bd 100644 --- a/src/pyclaw/classic/solver.py +++ b/src/pyclaw/classic/solver.py @@ -314,7 +314,7 @@ def step_hyperbolic(self,solution): # Find local value for dt/dx if state.index_capa>=0: - dtdx = self.dt / (grid.delta[0] * state.aux[state.index_capa,:]) + dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) else: dtdx += self.dt/grid.delta[0] From 688587035fe43d32627cb32f6d8007739affe8ee Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 13:39:55 -0400 Subject: [PATCH 07/47] Update advection_1d_nonunif.py fixed aux array --- examples/advection_1d/advection_1d_nonunif.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 9d7c83724..bcf7fbba9 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -25,7 +25,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, - time_integrator='SSP104', outdir='./_output_hb'): + time_integrator='SSP104', outdir='./_output'): if use_petsc: import clawpack.petclaw as pyclaw @@ -58,20 +58,16 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.aux_bc_lower[0] = pyclaw.BC.periodic solver.aux_bc_upper[0] = pyclaw.BC.periodic - x = pyclaw.Dimension(0.0,1.0,nx,name='x') + x = pyclaw.Dimension(0.0,1.0,nx,name='x',num_aux=1) domain = pyclaw.Domain(x) state = pyclaw.State(domain,1) state.problem_data['u'] = 1. # Advection velocity - nw = 50 # location of nouniformity - alpha = 0.1 # ratio of nonuniform small cell to standard cell state.index_capa = 0 - xpxc = nx * 1.0 / (nx-1) state.aux = np.zeros((1,nx)) - state.aux[0, :] = xpxc - state.aux[0, nw-1] = alpha * xpxc - state.aux[0, nw] = (1-alpha) * xpxc + state.aux[0, :int(nw)/2 + 1] = 0.9 + state.aux[0, int(nw)/2 + 1:] = 1.1 # Initial data xc = state.grid.x.centers @@ -83,7 +79,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi claw.solution = pyclaw.Solution(state,domain) claw.solver = solver - claw.tfinal =1.0 + claw.tfinal = 1.0 claw.outdir = outdir if outdir is None: claw.output_format = None From 0fa6e21c3b7c962f026c40224581393217a4c646 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 13:40:25 -0400 Subject: [PATCH 08/47] Update test_advection_nonunif.py fixed import --- examples/advection_1d/test_advection_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index f0e439a74..e8d586f9b 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -4,7 +4,7 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ - from . import advection_1d + from . import advection_1d_nonunif def verify_expected(expected): """ given an expected value, returns a verification function """ From 76272076cbe268ede63121419b1f1a70f3485eb9 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 13:58:47 -0400 Subject: [PATCH 09/47] Update advection_1d_nonunif.py fixed num_aux --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index bcf7fbba9..8b2319ae3 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -58,9 +58,9 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.aux_bc_lower[0] = pyclaw.BC.periodic solver.aux_bc_upper[0] = pyclaw.BC.periodic - x = pyclaw.Dimension(0.0,1.0,nx,name='x',num_aux=1) + x = pyclaw.Dimension(0.0,1.0,nx,name='x') domain = pyclaw.Domain(x) - state = pyclaw.State(domain,1) + state = pyclaw.State(domain,1, num_aux=1) state.problem_data['u'] = 1. # Advection velocity From bb20f1fca5d9826b09c5b22c0302bb1513a54db8 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 14:30:39 -0400 Subject: [PATCH 10/47] Update advection_1d_nonunif.py fix aux variable and xc to match nonuniform grid. --- examples/advection_1d/advection_1d_nonunif.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 8b2319ae3..435ef7b22 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -11,7 +11,8 @@ q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. - +The grid is nonuniform as shown by the aux variable, with first half of the grid being 0.9*dx and second half being 1.1*dx, +where dx is the computational grid size i.e. 1/nx The initial condition is a Gaussian and the boundary conditions are periodic. The final solution is identical to the initial data because the wave has crossed the domain exactly once. @@ -66,11 +67,15 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi state.index_capa = 0 state.aux = np.zeros((1,nx)) - state.aux[0, :int(nw)/2 + 1] = 0.9 - state.aux[0, int(nw)/2 + 1:] = 1.1 + state.aux[0, :int(nx)/2 + 1] = 0.9 + state.aux[0, int(nx)/2 + 1:] = 1.1 # Initial data xc = state.grid.x.centers + delta_xc = (1.0-0.0)/nx + xc[:int(nx/2)+1] -= 0.05*delta_xc + xc[int(nx/2)+1:] += 0.05*delta_xc + beta = 100; gamma = 0; x0 = 0.75 state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) From fba66c44f11aafac4b208efb3d72625a55348392 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 14:43:14 -0400 Subject: [PATCH 11/47] Update advection_1d_nonunif.py typo in aux variable --- examples/advection_1d/advection_1d_nonunif.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 435ef7b22..5f220a177 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -67,8 +67,8 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi state.index_capa = 0 state.aux = np.zeros((1,nx)) - state.aux[0, :int(nx)/2 + 1] = 0.9 - state.aux[0, int(nx)/2 + 1:] = 1.1 + state.aux[0, :int(nx/2) + 1] = 0.9 + state.aux[0, int(nx/2) + 1:] = 1.1 # Initial data xc = state.grid.x.centers @@ -94,9 +94,9 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi claw.setplot = setplot -# plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) +# plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) - # return claw + return claw def setplot(plotdata): """ From 5717fdef3d6e64fdaa8e64bb75675baa754d6183 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Tue, 8 Oct 2019 15:21:14 -0400 Subject: [PATCH 12/47] Update advection_1d_nonunif.py fix num_output_times to 10 to avoid outputting more steps. --- examples/advection_1d/advection_1d_nonunif.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 5f220a177..a7cf94834 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -86,6 +86,9 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi claw.tfinal = 1.0 claw.outdir = outdir + claw.num_output_times = 10 + claw.nstepout = 1 + if outdir is None: claw.output_format = None From 9360ca15bc9c526fdbd8e67236dce046bbad0b0b Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Thu, 10 Oct 2019 21:16:24 -0700 Subject: [PATCH 13/47] Update advection_1d_nonunif.py Added mapc2p --- examples/advection_1d/advection_1d_nonunif.py | 35 ++++++++++--------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index a7cf94834..c1f4b5e0c 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -22,7 +22,19 @@ from clawpack import riemann from clawpack.pyclaw.plot import plot - +def mapc2p_nonunif(xc): + """ + moves the computational centers--which are centers of uniform grid from [0,1] + to centers of physical grid, whose cell size is 0.9*dx_computational from [0, 0.45] + and 1.1*dx_computational from [0.45, 1] + """ + nx = 100 # number of cells + shift = np.ones(nx) + shift[:int(nx/2)+1] *= -0.05 * (1-0)/nx + shift[int(nx/2)+1:] *= 0.05 * (1-0)/nx + xp = xc + shift + + return xp def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, @@ -50,10 +62,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language - solver.order=1 - solver.limiters = None - solver.num_eqn=1 - solver.num_waves=1 + solver.bc_lower[0] = pyclaw.BC.periodic solver.bc_upper[0] = pyclaw.BC.periodic solver.aux_bc_lower[0] = pyclaw.BC.periodic @@ -71,12 +80,8 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi state.aux[0, int(nx/2) + 1:] = 1.1 # Initial data - xc = state.grid.x.centers - delta_xc = (1.0-0.0)/nx - xc[:int(nx/2)+1] -= 0.05*delta_xc - xc[int(nx/2)+1:] += 0.05*delta_xc - beta = 100; gamma = 0; x0 = 0.75 + xc = state.grid.x.centers state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) claw = pyclaw.Controller() @@ -92,11 +97,8 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi if outdir is None: claw.output_format = None - claw.run() - - claw.setplot = setplot - +# claw.run() # uncomment these two lines to also plot # plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) return claw @@ -106,7 +108,7 @@ def setplot(plotdata): Plot solution using VisClaw. """ plotdata.clearfigures() # clear any old figures,axes,items data - + plotdata.mapc2p = mapc2p_nonunif plotfigure = plotdata.new_plotfigure(name='q', figno=1) # Set up for axes in this figure: @@ -120,7 +122,8 @@ def setplot(plotdata): plotitem.plotstyle = '-o' plotitem.color = 'b' plotitem.kwargs = {'linewidth':2,'markersize':5} - + plotitem.MappedGrid = True + return plotdata if __name__=="__main__": From 429c96dd1b60ef343b5d0ccac9dfa182488bf5ed Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Thu, 10 Oct 2019 21:17:43 -0700 Subject: [PATCH 14/47] Update test_advection_nonunif.py test classic solver only for advection_1d_nonunif --- examples/advection_1d/test_advection_nonunif.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index e8d586f9b..e4642ca22 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -5,6 +5,7 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ from . import advection_1d_nonunif + from . import advection_1d # we test the classic solver only for advection_1d_nonunif def verify_expected(expected): """ given an expected value, returns a verification function """ @@ -25,7 +26,7 @@ def advection_verify(claw): from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(3.203924e-04), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(4.47796411e-4), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) From 4739bae1915113add333b11727825baae258a4d2 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Thu, 10 Oct 2019 21:19:41 -0700 Subject: [PATCH 15/47] Update advection_1d_nonunif.py --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index c1f4b5e0c..972df21b0 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -11,8 +11,8 @@ q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. -The grid is nonuniform as shown by the aux variable, with first half of the grid being 0.9*dx and second half being 1.1*dx, -where dx is the computational grid size i.e. 1/nx +The physical grid is nonuniform as set by mapc2p by a shifting of centers to left in between [0, 0.45] +and shifting of centers to right in between [0.45,1] . The initial condition is a Gaussian and the boundary conditions are periodic. The final solution is identical to the initial data because the wave has crossed the domain exactly once. From 1c4a3da8083ee44badcabbca64a3801a8f5672d1 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Thu, 10 Oct 2019 21:22:06 -0700 Subject: [PATCH 16/47] Update test_advection_nonunif.py test only classic solver for advection_1d_nonunif --- examples/advection_1d/test_advection_nonunif.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index e4642ca22..49b723f06 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -30,15 +30,15 @@ def advection_verify(claw): kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) - sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected(1.163605e-05), + sharp_tests_rk = gen_variants(advection_1d.setup, verify_expected(1.163605e-05), kernel_languages=('Python','Fortran'), solver_type='sharpclaw',time_integrator='SSP104', outdir=None) - sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected(1.500727e-05), + sharp_tests_lmm = gen_variants(advection_1d.setup, verify_expected(1.500727e-05), kernel_languages=('Python','Fortran'), solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) - weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(7.489618e-06), + weno_tests = gen_variants(advection_1d.setup, verify_expected(7.489618e-06), kernel_languages=('Fortran',), solver_type='sharpclaw', time_integrator='SSP104', weno_order=17, outdir=None) From 301333e7aaf28866b30729af0103760116621477 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 19:11:20 -0400 Subject: [PATCH 17/47] Delete advection_1d_nonunif.py --- examples/advection_1d/advection_1d_nonunif.py | 132 ------------------ 1 file changed, 132 deletions(-) delete mode 100644 examples/advection_1d/advection_1d_nonunif.py diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py deleted file mode 100644 index 972df21b0..000000000 --- a/examples/advection_1d/advection_1d_nonunif.py +++ /dev/null @@ -1,132 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -r""" -One-dimensional advection -========================= - -Solve the linear advection equation: - -.. math:: - q_t + u q_x = 0. - -Here q is the density of some conserved quantity and u is the velocity. -The physical grid is nonuniform as set by mapc2p by a shifting of centers to left in between [0, 0.45] -and shifting of centers to right in between [0.45,1] . -The initial condition is a Gaussian and the boundary conditions are periodic. -The final solution is identical to the initial data because the wave has -crossed the domain exactly once. -""" -from __future__ import absolute_import -import numpy as np -from clawpack import riemann -from clawpack.pyclaw.plot import plot - -def mapc2p_nonunif(xc): - """ - moves the computational centers--which are centers of uniform grid from [0,1] - to centers of physical grid, whose cell size is 0.9*dx_computational from [0, 0.45] - and 1.1*dx_computational from [0.45, 1] - """ - nx = 100 # number of cells - shift = np.ones(nx) - shift[:int(nx/2)+1] *= -0.05 * (1-0)/nx - shift[int(nx/2)+1:] *= 0.05 * (1-0)/nx - xp = xc + shift - - return xp - - -def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, - time_integrator='SSP104', outdir='./_output'): - - if use_petsc: - import clawpack.petclaw as pyclaw - else: - from clawpack import pyclaw - - if kernel_language == 'Fortran': - riemann_solver = riemann.advection_1D - elif kernel_language == 'Python': - riemann_solver = riemann.advection_1D_py.advection_1D - - if solver_type=='classic': - solver = pyclaw.ClawSolver1D(riemann_solver) - elif solver_type=='sharpclaw': - solver = pyclaw.SharpClawSolver1D(riemann_solver) - solver.weno_order = weno_order - solver.time_integrator = time_integrator - if time_integrator == 'SSPLMMk3': - solver.lmm_steps = 5 - solver.check_lmm_cond = True - else: raise Exception('Unrecognized value of solver_type.') - - solver.kernel_language = kernel_language - - solver.bc_lower[0] = pyclaw.BC.periodic - solver.bc_upper[0] = pyclaw.BC.periodic - solver.aux_bc_lower[0] = pyclaw.BC.periodic - solver.aux_bc_upper[0] = pyclaw.BC.periodic - - x = pyclaw.Dimension(0.0,1.0,nx,name='x') - domain = pyclaw.Domain(x) - state = pyclaw.State(domain,1, num_aux=1) - - state.problem_data['u'] = 1. # Advection velocity - - state.index_capa = 0 - state.aux = np.zeros((1,nx)) - state.aux[0, :int(nx/2) + 1] = 0.9 - state.aux[0, int(nx/2) + 1:] = 1.1 - - # Initial data - beta = 100; gamma = 0; x0 = 0.75 - xc = state.grid.x.centers - state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) - - claw = pyclaw.Controller() - claw.keep_copy = True - claw.solution = pyclaw.Solution(state,domain) - claw.solver = solver - - claw.tfinal = 1.0 - claw.outdir = outdir - claw.num_output_times = 10 - claw.nstepout = 1 - - if outdir is None: - claw.output_format = None - - claw.setplot = setplot -# claw.run() # uncomment these two lines to also plot -# plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) - - return claw - -def setplot(plotdata): - """ - Plot solution using VisClaw. - """ - plotdata.clearfigures() # clear any old figures,axes,items data - plotdata.mapc2p = mapc2p_nonunif - plotfigure = plotdata.new_plotfigure(name='q', figno=1) - - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.ylimits = [-.2,1.0] - plotaxes.title = 'q' - - # Set up for item on these axes: - plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - plotitem.plot_var = 0 - plotitem.plotstyle = '-o' - plotitem.color = 'b' - plotitem.kwargs = {'linewidth':2,'markersize':5} - plotitem.MappedGrid = True - - return plotdata - -if __name__=="__main__": - #from clawpack.pyclaw.util import run_app_from_main - #output = run_app_from_main(setup,setplot) - setup() From ccbb54adc3eea1c9acac4d06a03cd7191b9c5f8d Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 16:12:03 -0700 Subject: [PATCH 18/47] new nonuniform grid advection example transformation x**2 --- examples/advection_1d/advection_1d_nonunif.py | 138 ++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 examples/advection_1d/advection_1d_nonunif.py diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py new file mode 100644 index 000000000..1fb8cbecf --- /dev/null +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +# encoding: utf-8 + +r""" +One-dimensional advection +========================= + +Solve the linear advection equation: + +.. math:: + q_t + u q_x = 0. + +Here q is the density of some conserved quantity and u is the velocity. + +The initial condition is a Gaussian and the boundary conditions are periodic. +The final solution is identical to the initial data because the wave has +crossed the domain exactly once. +""" +from __future__ import absolute_import +import numpy as np +from clawpack import riemann +from clawpack.pyclaw.plot import plot + + + + +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, + time_integrator='SSP104', outdir='./_output_hb'): + + if use_petsc: + import clawpack.petclaw as pyclaw + else: + from clawpack import pyclaw + import clawpack.pyclaw.geometry as Geom + + if kernel_language == 'Fortran': + riemann_solver = riemann.advection_1D + elif kernel_language == 'Python': + riemann_solver = riemann.advection_1D_py.advection_1D + + if solver_type=='classic': + solver = pyclaw.ClawSolver1D(riemann_solver) + elif solver_type=='sharpclaw': + solver = pyclaw.SharpClawSolver1D(riemann_solver) + solver.weno_order = weno_order + solver.time_integrator = time_integrator + if time_integrator == 'SSPLMMk3': + solver.lmm_steps = 5 + solver.check_lmm_cond = True + else: raise Exception('Unrecognized value of solver_type.') + + solver.kernel_language = kernel_language + solver.order=1 + solver.limiters = None + solver.num_eqn=1 + solver.num_waves=1 + solver.bc_lower[0] = pyclaw.BC.periodic + solver.bc_upper[0] = pyclaw.BC.periodic + solver.aux_bc_lower[0] = pyclaw.BC.periodic + solver.aux_bc_upper[0] = pyclaw.BC.periodic + + x = Geom.Dimension(0.0,1.0,nx,name='x') + domain = Geom.Domain(x) + grid1d = Geom.Grid(x) + state = pyclaw.State(domain,1,num_aux=1) + state.problem_data['u'] = 1. # Advection velocity + + # mapping to nonunif grid + double = lambda xarr : np.array([x**2 for x in xarr]) + grid1d.mapc2p = double +# print(grid1d.p_centers) + state.aux = np.zeros((1,nx)) + for i in range(nx): + state.aux[0, i] = (2*i + 1)/nx # dxp_i/dxc_i + # dxp_i = (2i+1)/(nx^2) + # dxc_i = 1/nx + + # Initial data + xc = state.grid.x.centers + beta = 100; gamma = 0; x0 = 0.75 + state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) + + claw = pyclaw.Controller() + claw.keep_copy = True + claw.solution = pyclaw.Solution(state,domain) + claw.solver = solver + + claw.tfinal =1.0 + claw.outdir = outdir + claw.num_output_times = 10 + claw.nstepout = 1 + if outdir is None: + claw.output_format = None + +# claw.run() + + claw.setplot = setplot +# plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) + + return claw + +def mapc2p_nonunif(xc): + """computational coord center is given by xc_i = (2i+1)/(2nx) + physical coord center is given by xp_i = (i^2+i+1/2)/(nx^2) + to go from xc to xp, we isolate the i from xc_i and transform it according + to xp_i formula + """ + nx = 100 + I = (xc * 2 * nx - 1)/2 + xp = (1/nx)**2 * (np.square(I) + I + 1/2) + return xp + +def setplot(plotdata): + """ + Plot solution using VisClaw. + """ + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.mapc2p = mapc2p_nonunif + plotfigure = plotdata.new_plotfigure(name='q', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.ylimits = [-.2,1.0] + plotaxes.title = 'q' + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = '-o' + plotitem.color = 'b' + plotitem.kwargs = {'linewidth':2,'markersize':5} + + return plotdata + +if __name__=="__main__": + #from clawpack.pyclaw.util import run_app_from_main + #output = run_app_from_main(setup,setplot) + setup() From 358d649e8ea563ba90a69c48c1bc2223919ae951 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 19:14:11 -0400 Subject: [PATCH 19/47] Update advection_1d_nonunif.py add limiters for test file to run without error --- examples/advection_1d/advection_1d_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 1fb8cbecf..bf525bd49 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -51,7 +51,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.kernel_language = kernel_language solver.order=1 - solver.limiters = None + solver.limiters = pyclaw.limiters.tvd.minmod solver.num_eqn=1 solver.num_waves=1 solver.bc_lower[0] = pyclaw.BC.periodic From 4d381a3caf71393a5cb7cefbfd1e11322b92d105 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 19:16:06 -0400 Subject: [PATCH 20/47] Delete test_advection_nonunif.py --- .../advection_1d/test_advection_nonunif.py | 53 ------------------- 1 file changed, 53 deletions(-) delete mode 100644 examples/advection_1d/test_advection_nonunif.py diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py deleted file mode 100644 index 49b723f06..000000000 --- a/examples/advection_1d/test_advection_nonunif.py +++ /dev/null @@ -1,53 +0,0 @@ -from __future__ import absolute_import -def test_1d_advection(): - """test_1d_advection - - tests against expected classic, sharpclaw, and high-order weno results """ - - from . import advection_1d_nonunif - from . import advection_1d # we test the classic solver only for advection_1d_nonunif - - def verify_expected(expected): - """ given an expected value, returns a verification function """ - def advection_verify(claw): - from clawpack.pyclaw.util import check_diff - import numpy as np - - q0=claw.frames[0].state.get_q_global() - qfinal=claw.frames[claw.num_output_times].state.get_q_global() - - if q0 is not None and qfinal is not None: - dx=claw.solution.domain.grid.delta[0] - test = dx*np.linalg.norm(qfinal-q0,1) - return check_diff(expected, test, reltol=1e-4) - else: - return - return advection_verify - - from clawpack.pyclaw.util import gen_variants - - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(4.47796411e-4), - kernel_languages=('Python','Fortran'), - solver_type='classic', outdir=None) - - sharp_tests_rk = gen_variants(advection_1d.setup, verify_expected(1.163605e-05), - kernel_languages=('Python','Fortran'), - solver_type='sharpclaw',time_integrator='SSP104', outdir=None) - - sharp_tests_lmm = gen_variants(advection_1d.setup, verify_expected(1.500727e-05), - kernel_languages=('Python','Fortran'), - solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) - - weno_tests = gen_variants(advection_1d.setup, verify_expected(7.489618e-06), - kernel_languages=('Fortran',), solver_type='sharpclaw', - time_integrator='SSP104', weno_order=17, - outdir=None) - - from itertools import chain - for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests): - yield test - - -if __name__=="__main__": - import nose - nose.main() From e97de2570c853600edcfc007bbf20b703130c933 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 16:16:35 -0700 Subject: [PATCH 21/47] Add files via upload test file for nonuniform example: tests only classic solver for nonuniform example --- .../advection_1d/test_advection_nonunif.py | 76 +++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 examples/advection_1d/test_advection_nonunif.py diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py new file mode 100644 index 000000000..871f59b1d --- /dev/null +++ b/examples/advection_1d/test_advection_nonunif.py @@ -0,0 +1,76 @@ +from __future__ import absolute_import +def test_1d_advection(): + """test_1d_advection + + tests against expected classic, sharpclaw, and high-order weno results """ + + from . import advection_1d_nonunif + from . import advection_1d + from clawpack.pyclaw.util import check_diff + import numpy as np + + def verify_expected(expected): + """ given an expected value, returns a verification function """ + def advection_verify(claw): + from clawpack.pyclaw.util import check_diff + import numpy as np + + q0=claw.frames[0].state.get_q_global() + qfinal=claw.frames[claw.num_output_times].state.get_q_global() + + if q0 is not None and qfinal is not None: + dx=claw.solution.domain.grid.delta[0] + test = dx*np.linalg.norm(qfinal-q0,1) + print(test) + return check_diff(expected, test, reltol=1e-4) + else: + return + return advection_verify + + + def verify_expected_nu(expected): + # given an expected value, returns a verification function for nonuniform 1d advection example + + def advection_verify_nu(claw): + q0=claw.frames[0].state.get_q_global() + qfinal=claw.frames[claw.num_output_times].state.get_q_global() + if q0 is not None and qfinal is not None: + dx=claw.solution.domain.grid.delta[0] + nx = 100 + cap_arr = np.zeros(nx) + for i in range(nx): + cap_arr[i] = (2*i + 1)/nx + test = np.sum(np.multiply(dx*cap_arr,(qfinal-q0))) + return check_diff(expected, test, reltol=1e-4) + else: + return + return advection_verify_nu + + + from clawpack.pyclaw.util import gen_variants + + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-1.6905458739204524e-04), + kernel_languages=('Python','Fortran'), + solver_type='classic', outdir=None) + + sharp_tests_rk = gen_variants(advection_1d.setup, verify_expected(1.163605e-05), + kernel_languages=('Python','Fortran'), + solver_type='sharpclaw',time_integrator='SSP104', outdir=None) + + sharp_tests_lmm = gen_variants(advection_1d.setup, verify_expected(1.500727e-05), + kernel_languages=('Python','Fortran'), + solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) + + weno_tests = gen_variants(advection_1d.setup, verify_expected(7.489618e-06), + kernel_languages=('Fortran',), solver_type='sharpclaw', + time_integrator='SSP104', weno_order=17, + outdir=None) + + from itertools import chain + for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests): + yield test + + +if __name__=="__main__": + import nose + nose.main() From 4e79337cd601b8fa30ef223ef36b238a145ada67 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 19:18:02 -0400 Subject: [PATCH 22/47] Update advection_1d_nonunif.py added docstring --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index bf525bd49..963ad469e 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -5,13 +5,13 @@ One-dimensional advection ========================= -Solve the linear advection equation: +Solve the linear advection equation on a nonuniform grid: .. math:: q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. - +The nonuniform grid is accomplished by transformation of x**2 The initial condition is a Gaussian and the boundary conditions are periodic. The final solution is identical to the initial data because the wave has crossed the domain exactly once. From d0470df37cc6777bd1022e54bd3c82f2f57b410f Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:02:01 -0400 Subject: [PATCH 23/47] Update test_advection_nonunif.py fix test norm --- examples/advection_1d/test_advection_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 871f59b1d..9f0ed64fc 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -49,7 +49,7 @@ def advection_verify_nu(claw): from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-1.6905458739204524e-04), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(1.7746546860727922e-04), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) From de5f6f687264c577844a29aa1d959120b29fc84a Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:04:11 -0400 Subject: [PATCH 24/47] Update advection_1d_nonunif.py correct geometry setup --- examples/advection_1d/advection_1d_nonunif.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 963ad469e..1c4716bb4 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -31,7 +31,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi import clawpack.petclaw as pyclaw else: from clawpack import pyclaw - import clawpack.pyclaw.geometry as Geom + if kernel_language == 'Fortran': riemann_solver = riemann.advection_1D @@ -59,9 +59,9 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.aux_bc_lower[0] = pyclaw.BC.periodic solver.aux_bc_upper[0] = pyclaw.BC.periodic - x = Geom.Dimension(0.0,1.0,nx,name='x') - domain = Geom.Domain(x) - grid1d = Geom.Grid(x) + x = pyclaw.Dimension(0.0,1.0,nx,name='x') + domain = pyclaw.Domain(x) + grid1d = pyclaw.geometry.Grid(x) state = pyclaw.State(domain,1,num_aux=1) state.problem_data['u'] = 1. # Advection velocity From 4004e4a00611f8946d003d35b2407c946aca0fd1 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:18:46 -0400 Subject: [PATCH 25/47] Update test_advection_nonunif.py fix test norm --- examples/advection_1d/test_advection_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 9f0ed64fc..1712d9f6a 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -49,7 +49,7 @@ def advection_verify_nu(claw): from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(1.7746546860727922e-04), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-1.7746546860727922e-04), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) From ed3ad5917fb3e621ecbc685568ded782afca3a3f Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:20:25 -0400 Subject: [PATCH 26/47] Update advection_1d_nonunif.py fix geometry setup --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 1c4716bb4..5a5cc530e 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -31,7 +31,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi import clawpack.petclaw as pyclaw else: from clawpack import pyclaw - + from clawpack.pyclaw.geometry import Grid if kernel_language == 'Fortran': riemann_solver = riemann.advection_1D @@ -61,7 +61,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi x = pyclaw.Dimension(0.0,1.0,nx,name='x') domain = pyclaw.Domain(x) - grid1d = pyclaw.geometry.Grid(x) + grid1d = Grid(x) state = pyclaw.State(domain,1,num_aux=1) state.problem_data['u'] = 1. # Advection velocity From f92dae57f5aa69618fc279fd063d748499487c23 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:23:43 -0400 Subject: [PATCH 27/47] Update test_advection_nonunif.py --- examples/advection_1d/test_advection_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 1712d9f6a..a814197a9 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -49,7 +49,7 @@ def advection_verify_nu(claw): from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-1.7746546860727922e-04), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-0.00016905458739204524), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) From c1992cd39ceeb662ce01667fda07fac13ccdd819 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:42:12 -0400 Subject: [PATCH 28/47] Update advection_1d_nonunif.py --- examples/advection_1d/advection_1d_nonunif.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 5a5cc530e..963ad469e 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -31,7 +31,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi import clawpack.petclaw as pyclaw else: from clawpack import pyclaw - from clawpack.pyclaw.geometry import Grid + import clawpack.pyclaw.geometry as Geom if kernel_language == 'Fortran': riemann_solver = riemann.advection_1D @@ -59,9 +59,9 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.aux_bc_lower[0] = pyclaw.BC.periodic solver.aux_bc_upper[0] = pyclaw.BC.periodic - x = pyclaw.Dimension(0.0,1.0,nx,name='x') - domain = pyclaw.Domain(x) - grid1d = Grid(x) + x = Geom.Dimension(0.0,1.0,nx,name='x') + domain = Geom.Domain(x) + grid1d = Geom.Grid(x) state = pyclaw.State(domain,1,num_aux=1) state.problem_data['u'] = 1. # Advection velocity From b9bb602537af748c7e86386ed7e7cc377d0534f2 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 21:42:48 -0400 Subject: [PATCH 29/47] Update test_advection_nonunif.py --- examples/advection_1d/test_advection_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index a814197a9..871f59b1d 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -49,7 +49,7 @@ def advection_verify_nu(claw): from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-0.00016905458739204524), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-1.6905458739204524e-04), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) From 81a881b122a96958628a36ac763e77474a3ae911 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 22:16:42 -0400 Subject: [PATCH 30/47] Update advection_1d_nonunif.py fix petsc option for test file --- examples/advection_1d/advection_1d_nonunif.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 963ad469e..93320fd9c 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -27,11 +27,11 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output_hb'): - if use_petsc: - import clawpack.petclaw as pyclaw - else: - from clawpack import pyclaw - import clawpack.pyclaw.geometry as Geom + #if use_petsc: + import clawpack.petclaw as pyclaw + #else: + from clawpack import pyclaw + import clawpack.pyclaw.geometry as Geom if kernel_language == 'Fortran': riemann_solver = riemann.advection_1D From 05bd49093ce8da6ab2e2ae832d92ae4d3b4e3ba7 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 22:17:43 -0400 Subject: [PATCH 31/47] Update advection_1d_nonunif.py --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 93320fd9c..6a6b80912 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -27,8 +27,8 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output_hb'): - #if use_petsc: - import clawpack.petclaw as pyclaw + if use_petsc: + import clawpack.petclaw as pyclaw #else: from clawpack import pyclaw import clawpack.pyclaw.geometry as Geom From 2c098a619ff598c94837cc0f18d83c97489ecf88 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 22:48:28 -0400 Subject: [PATCH 32/47] Update test_advection_nonunif.py fix verification function to match python 2 --- examples/advection_1d/test_advection_nonunif.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 871f59b1d..7dbc68b88 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -30,22 +30,21 @@ def advection_verify(claw): def verify_expected_nu(expected): # given an expected value, returns a verification function for nonuniform 1d advection example - - def advection_verify_nu(claw): + def advection_verify_nu(claw): q0=claw.frames[0].state.get_q_global() - qfinal=claw.frames[claw.num_output_times].state.get_q_global() + qfinal=claw.frames[claw.num_output_times].state.get_q_glo$ if q0 is not None and qfinal is not None: dx=claw.solution.domain.grid.delta[0] nx = 100 - cap_arr = np.zeros(nx) - for i in range(nx): - cap_arr[i] = (2*i + 1)/nx + xc_edges = np.linspace(0.0,1.0,nx+1) + xp_edges = np.square(xc_edges) + cap_arr = np.diff(xp_edges)/np.diff(xc_edges) test = np.sum(np.multiply(dx*cap_arr,(qfinal-q0))) return check_diff(expected, test, reltol=1e-4) else: return return advection_verify_nu - + from clawpack.pyclaw.util import gen_variants From e97f317b3c15b48803f0ac8187095c9622a5a071 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 22:57:12 -0400 Subject: [PATCH 33/47] Update test_advection_nonunif.py typo --- examples/advection_1d/test_advection_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 7dbc68b88..2f8ea2d9a 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -32,7 +32,7 @@ def verify_expected_nu(expected): # given an expected value, returns a verification function for nonuniform 1d advection example def advection_verify_nu(claw): q0=claw.frames[0].state.get_q_global() - qfinal=claw.frames[claw.num_output_times].state.get_q_glo$ + qfinal=claw.frames[claw.num_output_times].state.get_q_global() if q0 is not None and qfinal is not None: dx=claw.solution.domain.grid.delta[0] nx = 100 From 5d7c309a7fa326179ef3a9bdcbf95fcfa244433d Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Sun, 13 Oct 2019 23:27:29 -0400 Subject: [PATCH 34/47] Update test_advection_nonunif.py indent error --- examples/advection_1d/test_advection_nonunif.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 2f8ea2d9a..47ad6b0ec 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -27,10 +27,9 @@ def advection_verify(claw): return return advection_verify - def verify_expected_nu(expected): # given an expected value, returns a verification function for nonuniform 1d advection example - def advection_verify_nu(claw): + def advection_verify_nu(claw): q0=claw.frames[0].state.get_q_global() qfinal=claw.frames[claw.num_output_times].state.get_q_global() if q0 is not None and qfinal is not None: From 9a0612628e58c03097f0bc5360360581e0c2cb61 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Thu, 17 Oct 2019 23:38:32 -0400 Subject: [PATCH 35/47] added nonuniform example files --- examples/advection_1d/advection_1d.py | 15 ++-- examples/advection_1d/advection_1d_nonunif.py | 71 +++++++++++-------- .../advection_1d/test_advection_nonunif.py | 38 ++++++++-- 3 files changed, 82 insertions(+), 42 deletions(-) diff --git a/examples/advection_1d/advection_1d.py b/examples/advection_1d/advection_1d.py index 10d661c2f..a86c9b444 100755 --- a/examples/advection_1d/advection_1d.py +++ b/examples/advection_1d/advection_1d.py @@ -25,7 +25,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, - time_integrator='SSP104', outdir='./_output_hb'): + time_integrator='SSP104', outdir='./_output'): if use_petsc: import clawpack.petclaw as pyclaw @@ -55,8 +55,6 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.num_waves=1 solver.bc_lower[0] = pyclaw.BC.periodic solver.bc_upper[0] = pyclaw.BC.periodic - solver.aux_bc_lower[0] = pyclaw.BC.periodic - solver.aux_bc_upper[0] = pyclaw.BC.periodic x = pyclaw.Dimension(0.0,1.0,nx,name='x') domain = pyclaw.Domain(x) @@ -79,13 +77,15 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi if outdir is None: claw.output_format = None - claw.run() + # claw.run() claw.setplot = setplot +# plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) - # return claw + + return claw def setplot(plotdata): """ @@ -109,7 +109,10 @@ def setplot(plotdata): return plotdata -#plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) +#plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) + + + if __name__=="__main__": #from clawpack.pyclaw.util import run_app_from_main diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 9d7c83724..73569336f 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -11,7 +11,7 @@ q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. - +Here we have a nonuniform grid, given by the transformation x**2 The initial condition is a Gaussian and the boundary conditions are periodic. The final solution is identical to the initial data because the wave has crossed the domain exactly once. @@ -21,16 +21,25 @@ from clawpack import riemann from clawpack.pyclaw.plot import plot +def mapc2p_nonunif(xc): + """computational coord center is given by xc_i = (2i+1)/(2nx) + physicachange change + """ + neg = -1*(xc < 0) + (xc > 0) + xp = xc**2 + xp = neg*xp + return xp -def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, - time_integrator='SSP104', outdir='./_output_hb'): +def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpclaw', weno_order=5, + time_integrator='SSP104', outdir='./_output'): if use_petsc: import clawpack.petclaw as pyclaw - else: - from clawpack import pyclaw + #else: + from clawpack import pyclaw + import clawpack.pyclaw.geometry if kernel_language == 'Fortran': riemann_solver = riemann.advection_1D @@ -49,8 +58,8 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language - solver.order=1 - solver.limiters = None + #solver.order = 1 + solver.limiters = pyclaw.tvd.minmod solver.num_eqn=1 solver.num_waves=1 solver.bc_lower[0] = pyclaw.BC.periodic @@ -58,55 +67,57 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi solver.aux_bc_lower[0] = pyclaw.BC.periodic solver.aux_bc_upper[0] = pyclaw.BC.periodic - x = pyclaw.Dimension(0.0,1.0,nx,name='x') + x = pyclaw.Dimension(-1.0,1.0,nx,name='x') domain = pyclaw.Domain(x) - state = pyclaw.State(domain,1) - + state = pyclaw.State(domain,1,num_aux=1) state.problem_data['u'] = 1. # Advection velocity - nw = 50 # location of nouniformity - alpha = 0.1 # ratio of nonuniform small cell to standard cell - state.index_capa = 0 - xpxc = nx * 1.0 / (nx-1) + xc = state.grid.x.centers + grid1d = state.grid +# print("pre",xc) + # mapping to nonunif grid + grid1d.mapc2p = mapc2p_nonunif + state.aux = np.zeros((1,nx)) - state.aux[0, :] = xpxc - state.aux[0, nw-1] = alpha * xpxc - state.aux[0, nw] = (1-alpha) * xpxc + state.aux[0,:] = np.diff(grid1d.p_nodes)/np.diff(state.grid.x.nodes) + + # for i in range(nx): + # state.aux[0, i] = (2*i + 1)/nx # dxp_i/dxc_i + # # dxp_i = (2i+1)/(nx^2) + # # dxc_i = 1/nx # Initial data - xc = state.grid.x.centers beta = 100; gamma = 0; x0 = 0.75 - state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) + state.q[0,:] = np.exp(-beta * (grid1d.p_centers-x0)**2) * np.cos(gamma * (grid1d.p_centers - x0)) claw = pyclaw.Controller() claw.keep_copy = True claw.solution = pyclaw.Solution(state,domain) claw.solver = solver - claw.tfinal =1.0 + claw.tfinal = 2.0 claw.outdir = outdir + claw.num_output_times = 10 + claw.nstepout = 1 if outdir is None: claw.output_format = None - claw.run() - - claw.setplot = setplot - # plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) - # return claw + return claw def setplot(plotdata): """ Plot solution using VisClaw. """ plotdata.clearfigures() # clear any old figures,axes,items data - + plotdata.mapc2p = mapc2p_nonunif plotfigure = plotdata.new_plotfigure(name='q', figno=1) # Set up for axes in this figure: plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [-1.0,1.0] plotaxes.ylimits = [-.2,1.0] plotaxes.title = 'q' @@ -115,11 +126,9 @@ def setplot(plotdata): plotitem.plot_var = 0 plotitem.plotstyle = '-o' plotitem.color = 'b' - plotitem.kwargs = {'linewidth':2,'markersize':5} - + plotitem.kwargs = {'linewidth':2,'markersize':1} return plotdata if __name__=="__main__": - #from clawpack.pyclaw.util import run_app_from_main - #output = run_app_from_main(setup,setplot) - setup() + from clawpack.pyclaw.util import run_app_from_main + output = run_app_from_main(setup,setplot) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index f0e439a74..f56d4127c 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -5,6 +5,7 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ from . import advection_1d + from . import advection_1d_nonunif def verify_expected(expected): """ given an expected value, returns a verification function """ @@ -23,22 +24,49 @@ def advection_verify(claw): return return advection_verify + def verify_expected_nonunif(expected): + """ given an expected value, returns a verification function for the nonuniform advection ex""" + from clawpack.pyclaw.util import check_diff + import numpy as np + + def mapc2p_nonunif(xc): + neg = -1*(xc < 0) + (xc > 0) + xp = xc**2 + xp = neg*xp + return xp + + def advection_nu_verify(claw): + q0=claw.frames[0].state.get_q_global() + qfinal=claw.frames[claw.num_output_times].state.get_q_global() + if q0 is not None and qfinal is not None: + dx=claw.solution.domain.grid.delta[0] + grid1d = claw.frames[0].state.grid + grid1d.mapc2p = mapc2p_nonunif + + aux = np.zeros((1,500)) + aux[0,:] = np.diff(grid1d.p_nodes)/np.diff(grid1d.x.nodes) + test = abs(np.sum(dx*aux[0,:]*(qfinal-q0))) + return check_diff(expected, test, reltol=1e-4) + else: + return + return advection_nu_verify + from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(3.203924e-04), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.156102048093389e-05), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) - sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected(1.163605e-05), + sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.2345173603527554e-07), kernel_languages=('Python','Fortran'), solver_type='sharpclaw',time_integrator='SSP104', outdir=None) - sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected(1.500727e-05), + sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.712477331998265e-07), kernel_languages=('Python','Fortran'), solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) - weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected(7.489618e-06), - kernel_languages=('Fortran',), solver_type='sharpclaw', + weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(4.014894137206369e-08), + kernel_languages=('Fortran',), solver_type='sharpclaw', time_integrator='SSP104', weno_order=17, outdir=None) From 6cfb39c5f5aaea0bf9ec638768507ce82ac59608 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Thu, 17 Oct 2019 23:39:29 -0400 Subject: [PATCH 36/47] Delete advection_1d_nonunif.py --- examples/advection_1d/advection_1d_nonunif.py | 138 ------------------ 1 file changed, 138 deletions(-) delete mode 100644 examples/advection_1d/advection_1d_nonunif.py diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py deleted file mode 100644 index 6a6b80912..000000000 --- a/examples/advection_1d/advection_1d_nonunif.py +++ /dev/null @@ -1,138 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -r""" -One-dimensional advection -========================= - -Solve the linear advection equation on a nonuniform grid: - -.. math:: - q_t + u q_x = 0. - -Here q is the density of some conserved quantity and u is the velocity. -The nonuniform grid is accomplished by transformation of x**2 -The initial condition is a Gaussian and the boundary conditions are periodic. -The final solution is identical to the initial data because the wave has -crossed the domain exactly once. -""" -from __future__ import absolute_import -import numpy as np -from clawpack import riemann -from clawpack.pyclaw.plot import plot - - - - -def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, - time_integrator='SSP104', outdir='./_output_hb'): - - if use_petsc: - import clawpack.petclaw as pyclaw - #else: - from clawpack import pyclaw - import clawpack.pyclaw.geometry as Geom - - if kernel_language == 'Fortran': - riemann_solver = riemann.advection_1D - elif kernel_language == 'Python': - riemann_solver = riemann.advection_1D_py.advection_1D - - if solver_type=='classic': - solver = pyclaw.ClawSolver1D(riemann_solver) - elif solver_type=='sharpclaw': - solver = pyclaw.SharpClawSolver1D(riemann_solver) - solver.weno_order = weno_order - solver.time_integrator = time_integrator - if time_integrator == 'SSPLMMk3': - solver.lmm_steps = 5 - solver.check_lmm_cond = True - else: raise Exception('Unrecognized value of solver_type.') - - solver.kernel_language = kernel_language - solver.order=1 - solver.limiters = pyclaw.limiters.tvd.minmod - solver.num_eqn=1 - solver.num_waves=1 - solver.bc_lower[0] = pyclaw.BC.periodic - solver.bc_upper[0] = pyclaw.BC.periodic - solver.aux_bc_lower[0] = pyclaw.BC.periodic - solver.aux_bc_upper[0] = pyclaw.BC.periodic - - x = Geom.Dimension(0.0,1.0,nx,name='x') - domain = Geom.Domain(x) - grid1d = Geom.Grid(x) - state = pyclaw.State(domain,1,num_aux=1) - state.problem_data['u'] = 1. # Advection velocity - - # mapping to nonunif grid - double = lambda xarr : np.array([x**2 for x in xarr]) - grid1d.mapc2p = double -# print(grid1d.p_centers) - state.aux = np.zeros((1,nx)) - for i in range(nx): - state.aux[0, i] = (2*i + 1)/nx # dxp_i/dxc_i - # dxp_i = (2i+1)/(nx^2) - # dxc_i = 1/nx - - # Initial data - xc = state.grid.x.centers - beta = 100; gamma = 0; x0 = 0.75 - state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0)) - - claw = pyclaw.Controller() - claw.keep_copy = True - claw.solution = pyclaw.Solution(state,domain) - claw.solver = solver - - claw.tfinal =1.0 - claw.outdir = outdir - claw.num_output_times = 10 - claw.nstepout = 1 - if outdir is None: - claw.output_format = None - -# claw.run() - - claw.setplot = setplot -# plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) - - return claw - -def mapc2p_nonunif(xc): - """computational coord center is given by xc_i = (2i+1)/(2nx) - physical coord center is given by xp_i = (i^2+i+1/2)/(nx^2) - to go from xc to xp, we isolate the i from xc_i and transform it according - to xp_i formula - """ - nx = 100 - I = (xc * 2 * nx - 1)/2 - xp = (1/nx)**2 * (np.square(I) + I + 1/2) - return xp - -def setplot(plotdata): - """ - Plot solution using VisClaw. - """ - plotdata.clearfigures() # clear any old figures,axes,items data - plotdata.mapc2p = mapc2p_nonunif - plotfigure = plotdata.new_plotfigure(name='q', figno=1) - - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.ylimits = [-.2,1.0] - plotaxes.title = 'q' - - # Set up for item on these axes: - plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - plotitem.plot_var = 0 - plotitem.plotstyle = '-o' - plotitem.color = 'b' - plotitem.kwargs = {'linewidth':2,'markersize':5} - - return plotdata - -if __name__=="__main__": - #from clawpack.pyclaw.util import run_app_from_main - #output = run_app_from_main(setup,setplot) - setup() From ac2d76611b350c69ce3afd75e587a5a01a2f23fd Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Thu, 17 Oct 2019 23:39:38 -0400 Subject: [PATCH 37/47] Delete test_advection_nonunif.py --- .../advection_1d/test_advection_nonunif.py | 74 ------------------- 1 file changed, 74 deletions(-) delete mode 100644 examples/advection_1d/test_advection_nonunif.py diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py deleted file mode 100644 index 47ad6b0ec..000000000 --- a/examples/advection_1d/test_advection_nonunif.py +++ /dev/null @@ -1,74 +0,0 @@ -from __future__ import absolute_import -def test_1d_advection(): - """test_1d_advection - - tests against expected classic, sharpclaw, and high-order weno results """ - - from . import advection_1d_nonunif - from . import advection_1d - from clawpack.pyclaw.util import check_diff - import numpy as np - - def verify_expected(expected): - """ given an expected value, returns a verification function """ - def advection_verify(claw): - from clawpack.pyclaw.util import check_diff - import numpy as np - - q0=claw.frames[0].state.get_q_global() - qfinal=claw.frames[claw.num_output_times].state.get_q_global() - - if q0 is not None and qfinal is not None: - dx=claw.solution.domain.grid.delta[0] - test = dx*np.linalg.norm(qfinal-q0,1) - print(test) - return check_diff(expected, test, reltol=1e-4) - else: - return - return advection_verify - - def verify_expected_nu(expected): - # given an expected value, returns a verification function for nonuniform 1d advection example - def advection_verify_nu(claw): - q0=claw.frames[0].state.get_q_global() - qfinal=claw.frames[claw.num_output_times].state.get_q_global() - if q0 is not None and qfinal is not None: - dx=claw.solution.domain.grid.delta[0] - nx = 100 - xc_edges = np.linspace(0.0,1.0,nx+1) - xp_edges = np.square(xc_edges) - cap_arr = np.diff(xp_edges)/np.diff(xc_edges) - test = np.sum(np.multiply(dx*cap_arr,(qfinal-q0))) - return check_diff(expected, test, reltol=1e-4) - else: - return - return advection_verify_nu - - - from clawpack.pyclaw.util import gen_variants - - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nu(-1.6905458739204524e-04), - kernel_languages=('Python','Fortran'), - solver_type='classic', outdir=None) - - sharp_tests_rk = gen_variants(advection_1d.setup, verify_expected(1.163605e-05), - kernel_languages=('Python','Fortran'), - solver_type='sharpclaw',time_integrator='SSP104', outdir=None) - - sharp_tests_lmm = gen_variants(advection_1d.setup, verify_expected(1.500727e-05), - kernel_languages=('Python','Fortran'), - solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) - - weno_tests = gen_variants(advection_1d.setup, verify_expected(7.489618e-06), - kernel_languages=('Fortran',), solver_type='sharpclaw', - time_integrator='SSP104', weno_order=17, - outdir=None) - - from itertools import chain - for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests): - yield test - - -if __name__=="__main__": - import nose - nose.main() From b39106b7aba827451cc30eb539d3dc21d4ff17bc Mon Sep 17 00:00:00 2001 From: cr2940 Date: Thu, 17 Oct 2019 23:49:10 -0400 Subject: [PATCH 38/47] fixed up advection example files --- examples/advection_1d/advection_1d.py | 17 ++------------ examples/advection_1d/advection_1d_nonunif.py | 22 ++++++------------- 2 files changed, 9 insertions(+), 30 deletions(-) diff --git a/examples/advection_1d/advection_1d.py b/examples/advection_1d/advection_1d.py index a86c9b444..89ab557a1 100755 --- a/examples/advection_1d/advection_1d.py +++ b/examples/advection_1d/advection_1d.py @@ -22,8 +22,6 @@ from clawpack.pyclaw.plot import plot - - def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output'): @@ -77,14 +75,8 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi if outdir is None: claw.output_format = None - # claw.run() - - claw.setplot = setplot -# plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) - - return claw def setplot(plotdata): @@ -109,12 +101,7 @@ def setplot(plotdata): return plotdata -#plot(setplot=setplot,outdir='./_output',plotdir='./plots',iplot=False, htmlplot=True) - - - if __name__=="__main__": - #from clawpack.pyclaw.util import run_app_from_main - #output = run_app_from_main(setup,setplot) - setup() + from clawpack.pyclaw.util import run_app_from_main + output = run_app_from_main(setup,setplot) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 73569336f..9f1b0dfd2 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -11,19 +11,19 @@ q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. -Here we have a nonuniform grid, given by the transformation x**2 +Here we have a nonuniform grid, given by the transformation x**2 on grid [-1,1] The initial condition is a Gaussian and the boundary conditions are periodic. The final solution is identical to the initial data because the wave has -crossed the domain exactly once. +crossed the domain exactly once, which takes computational time 2, because the speed is 1 and grid length 2. """ from __future__ import absolute_import import numpy as np from clawpack import riemann -from clawpack.pyclaw.plot import plot def mapc2p_nonunif(xc): - """computational coord center is given by xc_i = (2i+1)/(2nx) - physicachange change + """This function takes the interval [-1,1] and squares the computational coordinate + while keeping the negative coordinates to be squared and yet retain their negativity + This serves as an example of transformation to a nonuniform grid. """ neg = -1*(xc < 0) + (xc > 0) xp = xc**2 @@ -37,7 +37,6 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc if use_petsc: import clawpack.petclaw as pyclaw - #else: from clawpack import pyclaw import clawpack.pyclaw.geometry @@ -58,7 +57,7 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language - #solver.order = 1 + #solver.order = 1 # in order to make test file pass without complaint solver.limiters = pyclaw.tvd.minmod solver.num_eqn=1 solver.num_waves=1 @@ -74,18 +73,12 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc xc = state.grid.x.centers grid1d = state.grid -# print("pre",xc) + # mapping to nonunif grid grid1d.mapc2p = mapc2p_nonunif - state.aux = np.zeros((1,nx)) state.aux[0,:] = np.diff(grid1d.p_nodes)/np.diff(state.grid.x.nodes) - # for i in range(nx): - # state.aux[0, i] = (2*i + 1)/nx # dxp_i/dxc_i - # # dxp_i = (2i+1)/(nx^2) - # # dxc_i = 1/nx - # Initial data beta = 100; gamma = 0; x0 = 0.75 state.q[0,:] = np.exp(-beta * (grid1d.p_centers-x0)**2) * np.cos(gamma * (grid1d.p_centers - x0)) @@ -103,7 +96,6 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc claw.output_format = None claw.setplot = setplot -# plot(setplot=setplot,outdir='./_output_hb',plotdir='./plots_hb',iplot=False, htmlplot=True) return claw From 609ed9b2ac9f009f9b642fbfd041edf34e694d80 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Thu, 17 Oct 2019 23:52:36 -0400 Subject: [PATCH 39/47] fixed verification function --- examples/advection_1d/test_advection_nonunif.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index f56d4127c..e33280b2a 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -42,8 +42,8 @@ def advection_nu_verify(claw): dx=claw.solution.domain.grid.delta[0] grid1d = claw.frames[0].state.grid grid1d.mapc2p = mapc2p_nonunif - - aux = np.zeros((1,500)) + nx = 500 + aux = np.zeros((1,nx)) aux[0,:] = np.diff(grid1d.p_nodes)/np.diff(grid1d.x.nodes) test = abs(np.sum(dx*aux[0,:]*(qfinal-q0))) return check_diff(expected, test, reltol=1e-4) From 04d3842171ff0a8f5434f077fe848dd3e22cfb43 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Fri, 18 Oct 2019 22:51:15 -0400 Subject: [PATCH 40/47] fixed up nonuniform advection example --- examples/advection_1d/__init__.py | 0 examples/advection_1d/advection_1d_nonunif.py | 28 +- src/pyclaw/classic/solverjli.py | 866 ++++++++++++++++++ src/pyclaw/sharpclaw/solver.py | 95 +- 4 files changed, 928 insertions(+), 61 deletions(-) delete mode 100644 examples/advection_1d/__init__.py create mode 100644 src/pyclaw/classic/solverjli.py diff --git a/examples/advection_1d/__init__.py b/examples/advection_1d/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 9f1b0dfd2..89567a0f8 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -11,18 +11,19 @@ q_t + u q_x = 0. Here q is the density of some conserved quantity and u is the velocity. -Here we have a nonuniform grid, given by the transformation x**2 on grid [-1,1] -The initial condition is a Gaussian and the boundary conditions are periodic. +Here we have a nonuniform grid, given by the transformation x**2 on grid [-0.5,0.5] to [-0.25,0.25]. +The initial condition is a Gaussian centered at 0 and the boundary conditions are periodic. The final solution is identical to the initial data because the wave has -crossed the domain exactly once, which takes computational time 2, because the speed is 1 and grid length 2. +crossed the domain exactly once, which takes computational time 0.5, because the speed is 1 and grid length 0.5. """ from __future__ import absolute_import import numpy as np from clawpack import riemann def mapc2p_nonunif(xc): - """This function takes the interval [-1,1] and squares the computational coordinate + """This function takes the interval [-xL,xR] and squares the computational coordinate while keeping the negative coordinates to be squared and yet retain their negativity + to the physical coordinate [-xL**2, xR**2] This serves as an example of transformation to a nonuniform grid. """ neg = -1*(xc < 0) + (xc > 0) @@ -32,8 +33,8 @@ def mapc2p_nonunif(xc): -def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpclaw', weno_order=5, - time_integrator='SSP104', outdir='./_output'): +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='sharpclaw', weno_order=5, + time_integrator='SSPLMMk3', outdir='./_output'): if use_petsc: import clawpack.petclaw as pyclaw @@ -57,7 +58,7 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language - #solver.order = 1 # in order to make test file pass without complaint + #solver.order = 2 # in order to make test file pass without complaint solver.limiters = pyclaw.tvd.minmod solver.num_eqn=1 solver.num_waves=1 @@ -66,21 +67,22 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc solver.aux_bc_lower[0] = pyclaw.BC.periodic solver.aux_bc_upper[0] = pyclaw.BC.periodic - x = pyclaw.Dimension(-1.0,1.0,nx,name='x') + x = pyclaw.Dimension(-0.5,0.5,nx,name='x') domain = pyclaw.Domain(x) state = pyclaw.State(domain,1,num_aux=1) state.problem_data['u'] = 1. # Advection velocity + state.index_capa = 0 xc = state.grid.x.centers grid1d = state.grid # mapping to nonunif grid grid1d.mapc2p = mapc2p_nonunif - state.aux = np.zeros((1,nx)) + state.aux = np.zeros((1,nx)) # capacity array dx_p/dx_c state.aux[0,:] = np.diff(grid1d.p_nodes)/np.diff(state.grid.x.nodes) # Initial data - beta = 100; gamma = 0; x0 = 0.75 + beta = 100; gamma = 0; x0 = 0.0 state.q[0,:] = np.exp(-beta * (grid1d.p_centers-x0)**2) * np.cos(gamma * (grid1d.p_centers - x0)) claw = pyclaw.Controller() @@ -88,7 +90,7 @@ def setup(nx=500, kernel_language='Python', use_petsc=False, solver_type='sharpc claw.solution = pyclaw.Solution(state,domain) claw.solver = solver - claw.tfinal = 2.0 + claw.tfinal = 0.5 # one cycle claw.outdir = outdir claw.num_output_times = 10 claw.nstepout = 1 @@ -109,7 +111,7 @@ def setplot(plotdata): # Set up for axes in this figure: plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = [-1.0,1.0] + plotaxes.xlimits = [-0.25,0.25] plotaxes.ylimits = [-.2,1.0] plotaxes.title = 'q' @@ -118,7 +120,7 @@ def setplot(plotdata): plotitem.plot_var = 0 plotitem.plotstyle = '-o' plotitem.color = 'b' - plotitem.kwargs = {'linewidth':2,'markersize':1} + plotitem.kwargs = {'linewidth':2,'markersize':5} return plotdata if __name__=="__main__": diff --git a/src/pyclaw/classic/solverjli.py b/src/pyclaw/classic/solverjli.py new file mode 100644 index 000000000..7dfd80c57 --- /dev/null +++ b/src/pyclaw/classic/solverjli.py @@ -0,0 +1,866 @@ +r""" +Module containing the classic Clawpack solvers. + +This module contains the pure and wrapped classic clawpack solvers. All +clawpack solvers inherit from the :class:`ClawSolver` superclass which in turn +inherits from the :class:`~pyclaw.solver.Solver` superclass. These +are both pure virtual classes; the only solver classes that should be instantiated +are the dimension-specific ones, :class:`ClawSolver1D` and :class:`ClawSolver2D`. +""" + +from clawpack.pyclaw.util import add_parent_doc +from clawpack.pyclaw.solver import Solver +from clawpack.pyclaw.limiters import tvd + +# ============================================================================ +# Generic Clawpack solver class +# ============================================================================ +class ClawSolver(Solver): + r""" + Generic classic Clawpack solver + + All Clawpack solvers inherit from this base class. + + .. attribute:: mthlim + + Limiter(s) to be used. Specified either as one value or a list. + If one value, the specified limiter is used for all wave families. + If a list, the specified values indicate which limiter to apply to + each wave family. Take a look at pyclaw.limiters.tvd for an enumeration. + ``Default = limiters.tvd.minmod`` + + .. attribute:: order + + Order of the solver, either 1 for first order (i.e., Godunov's method) + or 2 for second order (Lax-Wendroff-LeVeque). + ``Default = 2`` + + .. attribute:: source_split + + Which source splitting method to use: 1 for first + order Godunov splitting and 2 for second order Strang splitting. + ``Default = 1`` + + .. attribute:: fwave + + Whether to split the flux jump (rather than the jump in Q) into waves; + requires that the Riemann solver performs the splitting. + ``Default = False`` + + .. attribute:: step_source + + Handle for function that evaluates the source term. + The required signature for this function is: + + def step_source(solver,state,dt) + + .. attribute:: kernel_language + + Specifies whether to use wrapped Fortran routines ('Fortran') + or pure Python ('Python'). ``Default = 'Fortran'``. + + .. attribute:: verbosity + + The level of detail of logged messages from the Fortran solver. + ``Default = 0``. + + """ + + # ========== Generic Init Routine ======================================== + def __init__(self,riemann_solver=None,claw_package=None): + r""" + See :class:`ClawSolver` for full documentation. + + Output: + - (:class:`ClawSolver`) - Initialized clawpack solver + """ + self.num_ghost = 2 + self.limiters = tvd.vanleer + self.order = 2 + self.source_split = 1 + self.fwave = False + self.step_source = None + self.kernel_language = 'Fortran' + self.verbosity = 0 + self.cfl_max = 1.0 + self.cfl_desired = 0.9 + self._mthlim = self.limiters + self._method = None + self.dt_old = None + + # Call general initialization function + super(ClawSolver,self).__init__(riemann_solver,claw_package) + + # ========== Time stepping routines ====================================== + def step(self,solution,take_one_step,tstart,tend): + r""" + Evolve solution one time step + + The elements of the algorithm for taking one step are: + + 1. Pick a step size as specified by the base solver attribute :func:`get_dt` + + 2. A half step on the source term :func:`step_source` if Strang splitting is + being used (:attr:`source_split` = 2) + + 3. A step on the homogeneous problem :math:`q_t + f(q)_x = 0` is taken + + 4. A second half step or a full step is taken on the source term + :func:`step_source` depending on whether Strang splitting was used + (:attr:`source_split` = 2) or Godunov splitting + (:attr:`source_split` = 1) + + This routine is called from the method evolve_to_time defined in the + pyclaw.solver.Solver superclass. + + :Input: + - *solution* - (:class:`~pyclaw.solution.Solution`) solution to be evolved + + :Output: + - (bool) - True if full step succeeded, False otherwise + """ + self.get_dt(solution.t,tstart,tend,take_one_step) + self.cfl.set_global_max(0.) + + if self.source_split == 2 and self.step_source is not None: + self.step_source(self,solution.states[0],self.dt/2.0) + + self.step_hyperbolic(solution) + + # Check here if the CFL condition is satisfied. + # If not, return # immediately to evolve_to_time and let it deal with + # picking a new step size (dt). + if self.cfl.get_cached_max() >= self.cfl_max: + return False + + if self.step_source is not None: + # Strang splitting + if self.source_split == 2: + self.step_source(self,solution.states[0],self.dt/2.0) + + # Godunov Splitting + if self.source_split == 1: + self.step_source(self,solution.states[0],self.dt) + + return True + + def _check_cfl_settings(self): + pass + + def _allocate_workspace(self,solution): + pass + + def step_hyperbolic(self,solution): + r""" + Take one homogeneous step on the solution. + + This is a dummy routine and must be overridden. + """ + raise Exception("Dummy routine, please override!") + + def _set_mthlim(self): + r""" + Convenience routine to convert users limiter specification to + the format understood by the Fortran code (i.e., a list of length num_waves). + """ + #self._mthlim = self.limiters + #if not isinstance(self.limiters,list): self._mthlim=[self._mthlim] + # if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves + # if len(self._mthlim)!=self.num_waves: + # raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves') + + def _set_method(self,state): + r""" + Set values of the solver._method array required by the Fortran code. + These are algorithmic parameters. + """ + import numpy as np + #We ought to put method and many other things in a Fortran + #module and set the fortran variables directly here. + self._method =np.empty(7, dtype=int,order='F') + self._method[0] = self.dt_variable + self._method[1] = self.order + if self.num_dim==1: + self._method[2] = 0 # Not used in 1D + elif self.dimensional_split: + self._method[2] = -1 # First-order dimensional splitting + else: + self._method[2] = self.transverse_waves + self._method[3] = self.verbosity + self._method[4] = 0 # Not used for PyClaw (would be self.source_split) + self._method[5] = state.index_capa + 1 + self._method[6] = state.num_aux + + def setup(self,solution): + r""" + Perform essential solver setup. This routine must be called before + solver.step() may be called. + """ + # This is a hack to deal with the fact that petsc4py + # doesn't allow us to change the stencil_width (num_ghost) + solution.state.set_num_ghost(self.num_ghost) + # End hack + + self._check_cfl_settings() + + self._set_mthlim() + if(self.kernel_language == 'Fortran'): + if self.fmod is None: + so_name = 'clawpack.pyclaw.classic.classic'+str(self.num_dim) + self.fmod = __import__(so_name,fromlist=['clawpack.pyclaw.classic']) + self._set_fortran_parameters(solution) + self._allocate_workspace(solution) + elif self.num_dim>1: + raise Exception('Only Fortran kernels are supported in multi-D.') + + self._allocate_bc_arrays(solution.states[0]) + + super(ClawSolver,self).setup(solution) + + + def _set_fortran_parameters(self,solution): + r""" + Pack parameters into format recognized by Clawpack (Fortran) code. + + Sets the solver._method array and the cparam common block for the Riemann solver. + """ + self._set_method(solution.state) + # The reload here is necessary because otherwise the common block + # cparam in the Riemann solver doesn't get flushed between running + # different tests in a single Python session. + reload(self.fmod) + solution.state.set_cparam(self.fmod) + solution.state.set_cparam(self.rp) + + def __del__(self): + r""" + Delete Fortran objects, which otherwise tend to persist in Python sessions. + """ + if(self.kernel_language == 'Fortran'): + del self.fmod + + super(ClawSolver,self).__del__() + + +# ============================================================================ +# ClawPack 1d Solver Class +# ============================================================================ +class ClawSolver1D(ClawSolver): + r""" + Clawpack evolution routine in 1D + + This class represents the 1d clawpack solver on a single grid. Note that + there are routines here for interfacing with the fortran time stepping + routines and the Python time stepping routines. The ones used are + dependent on the argument given to the initialization of the solver + (defaults to python). + + """ + + __doc__ += add_parent_doc(ClawSolver) + + def __init__(self, riemann_solver=None, claw_package=None): + r""" + Create 1d Clawpack solver + + Output: + - (:class:`ClawSolver1D`) - Initialized 1d clawpack solver + + See :class:`ClawSolver1D` for more info. + """ + self.num_dim = 1 + self.reflect_index = [1] + + super(ClawSolver1D,self).__init__(riemann_solver, claw_package) + + + # ========== Homogeneous Step ===================================== + def step_hyperbolic(self,solution): + r""" + Take one time step on the homogeneous hyperbolic system. + + :Input: + - *solution* - (:class:`~pyclaw.solution.Solution`) Solution that + will be evolved + """ + import numpy as np + + state = solution.states[0] + grid = state.grid + + + self._apply_bcs(state) + + num_eqn,num_ghost = state.num_eqn,self.num_ghost + + if(self.kernel_language == 'Fortran'): + mx = grid.num_cells[0] + dx,dt = grid.delta[0],self.dt + dtdx = np.zeros( (mx+2*num_ghost) ) + dt/dx + rp1 = self.rp.rp1._cpointer + + self.qbc,cfl = self.fmod.step1(num_ghost,mx,self.qbc,self.auxbc,dx,dt,self._method,self._mthlim,self.fwave,rp1) + + elif(self.kernel_language == 'Python'): + q = self.qbc + aux = self.auxbc + # Limiter to use in the pth family + limiter = np.array(self._mthlim,ndmin=1) + + dtdx = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) + + # Find local value for dt/dx + if 'method' not in state.problem_data: + if state.index_capa>=0: + dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) + else: + dtdx += self.dt / grid.delta[0] + elif state.problem_data['method'] == 'h_box': + xpxc = state.problem_data['xpxc'] + dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) + dtdx_hbox = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) + dtdx_hbox += self.dt / (grid.delta[0] * xpxc) + elif state.problem_data['method'] == 'h_box_wave': + xpxc = state.problem_data['xpxc'] + dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) + dtdx_hbox = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) + dtdx_hbox += self.dt / (grid.delta[0] * xpxc) + + + + # Solve Riemann problem at each interface + q_l=q[:,:-1].copy() + q_r=q[:,1:].copy() + if state.aux is not None: + aux_l=aux[:,:-1].copy() + aux_r=aux[:,1:].copy() + else: + aux_l = None + aux_r = None + + if 'method' not in state.problem_data: + # normal case + wave,s,amdq,apdq = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) + + + # Update loop limits, these are the limits for the Riemann solver + # locations, which then update a grid cell value + # We include the Riemann problem just outside of the grid so we can + # do proper limiting at the grid edges + # LL | | UL + # | LL | | | | ... | | | UL | | + # | | + + LL = self.num_ghost - 1 + UL = self.num_ghost + grid.num_cells[0] + 1 + + # Update q for Godunov update + for m in range(num_eqn): + q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] + q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] + + elif state.problem_data['method'] == 'h_box': + # # add corrections + wave,s,amdq,apdq,f_corr_l,f_corr_r = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) + LL = self.num_ghost - 1 + UL = self.num_ghost + grid.num_cells[0] + 1 + + # Update q for Godunov update + for m in range(num_eqn): + q[m,LL:UL] -= dtdx[LL:UL]*(apdq[m,LL-1:UL-1] - f_corr_r[m,LL-1:UL-1]) + q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*(amdq[m,LL-1:UL-1] + f_corr_l[m,LL-1:UL-1]) + + elif state.problem_data['method'] == 'h_box_wave': + # # add corrections + state.problem_data['arrival_state'] = False + + + wave,s,amdq,apdq,q_hbox_initial,aux_hbox = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) + LL = self.num_ghost - 1 + UL = self.num_ghost + grid.num_cells[0] + 1 + + # Update q for Godunov update + iw = state.problem_data['wall_position'] + self.num_ghost - 1 + q_last = q[:,iw:iw+2].copy() + + for m in range(num_eqn): + q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] + q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] + + # check the arrivals + q[:,iw:iw+2] = q_last[:,:] # reset the wall cells + dt = self.dt + num_waves = self.num_waves + dx = grid.delta[0] * xpxc + alpha = state.problem_data['fraction'] + arrival_times = np.array([0.0]) + for mw in range(num_waves): + if (s[mw,iw-1] > 0 and (s[mw,iw-1] * dt > alpha * dx)): + arrival_times = np.append(arrival_times, alpha * dx / s[mw,iw-1]) + if (s[mw,iw+1] < 0 and ( (-s[mw,iw+1]) * dt > (1 - alpha) * dx ) ): + arrival_times = np.append(arrival_times, -(1 - alpha) * dx / s[mw,iw+1]) + arrival_times.sort() + + n_arrival_times = len(arrival_times) + if n_arrival_times == 1 : + state.problem_data['arrival_state'] = False + else: + state.problem_data['arrival_state'] = True + s_cells = np.zeros((num_waves, 3, n_arrival_times)) + s_cells[:,:,0] = s[:, iw-1:iw+2].copy() + wave_cells = np.zeros((num_eqn, num_waves, 3, n_arrival_times)) + wave_cells[:,:,:,0] = wave[:,:,iw-1:iw+2].copy() + + + if state.problem_data['arrival_state'] == False: + q[:,iw] -= dt/(alpha * dx) * apdq[:,iw-1] + q[:,iw+1] -= dt/((1 - alpha)*dx) * amdq[:,iw+1] + for mw in range(num_waves): + if (s[mw,iw] < 0): + q[:,iw-1] -= dt/dx * ( max(0, -s[mw,iw] * dt - alpha * dx) / (-s[mw,iw] * dt) * wave[:,mw,iw] ) + q[:,iw] -= dt/(alpha * dx) * (min(-s[mw,iw] * dt, alpha * dx) / (-s[mw,iw] * dt) * wave[:,mw,iw] ) + elif (s[mw,iw] > 0): + q[:,iw+1] -= dt/((1 - alpha)*dx) * (min(s[mw,iw] * dt, (1 - alpha) * dx) / (s[mw,iw] * dt) * wave[:,mw,iw] ) + q[:,iw+2] -= dt/dx * ( (max(0, s[mw,iw] * dt - (1 - alpha) * dx) / s[mw,iw] * dt) * wave[:,mw,iw] ) + + if state.problem_data['arrival_state'] == True: + ## update q_hbox + for i in range(1, n_arrival_times): + q_hbox = q_hbox_initial.copy() + + for mw in range(num_waves): + if s[mw,iw-2] > 0: + q_hbox[:,0] -= arrival_times[i] / dx * (max(0, s[mw,iw-2] * arrival_times[i] - alpha * dx) / (s[mw,iw-2] * arrival_times[i]) * wave[:,mw,iw-2]) + + if s[mw, iw-1] < 0: + q_hbox[:,0] -= arrival_times[i] / dx * (min(-s[mw,iw-1] * arrival_times[i], (1 - alpha) * dx) / (-s[mw,iw-1] * arrival_times[i]) * wave[:,mw,iw-1]) + + if s_cells[mw,0,i] > 0: + for j in range(i): + q_hbox[:,0] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,0,j]) + + if s_cells[mw,0,i] * arrival_times[i] > alpha * dx: # - 1e-14: + # check the arrival wave + wave_cells[:,mw,0,i] = 0.0 + + if s_cells[mw,1,i] < 0: + for j in range(i): + q_hbox[:,0] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,1,j]) + + if s_cells[mw,1,i] > 0: + for j in range(i): + q_hbox[:,1] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,1,j]) + + if s_cells[mw,2,i] < 0: + for j in range(i): + q_hbox[:,1] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,2,j]) + + if (-s_cells[mw,2,i] * arrival_times[i]) > (1 - alpha) * dx: # - 1e-14: + # check the arrival wave + wave_cells[:,mw,2,i] = 0.0 + + if s[mw,iw+1] > 0: + q_hbox[:,1] -= arrival_times[i] / dx * (min(s[mw,iw+1] * arrival_times[i], alpha * dx) / (-s[mw,iw+1] * arrival_times[i]) * wave[:,mw,iw+1]) + + if s[mw,iw+2] < 0: + q_hbox[:,1] -= arrival_times[i] / dx * (max(0, -s[mw,iw+2] * arrival_times[i] - (1 - alpha) * dx) / (-s[mw,iw+2] * arrival_times[i]) * wave[:,mw,iw+2]) + + + + + wave_cells[:,:,1,i],s_cells[:,1,i],amdq_arr,apdq_arr = self.rp(q_hbox[:,0],q_hbox[:,1],aux_hbox[:,0],aux_hbox[:,1],state.problem_data) +#update wave cells :,:,0,i and others by doing middle value update every arrival step. make qhbox_middle +#add code in shallow_fwave_hbox_dry_1d with fw[iw-1] = middle and q_iw-1 and fw[iw+1] = middle and q_iw+2 + ## update q[iw-1], q[iw], q[iw+1] and q[iw+2] + arrival_times = np.append(arrival_times, dt) + n_arrival_times = len(arrival_times) + + for mw in range(num_waves): + for i in range(n_arrival_times-1): + if s_cells[mw,0,i] > 0: + q[:,iw] -= (arrival_times[i+1] - arrival_times[i]) / (alpha * dx) * (wave_cells[:,mw,0,i]) + if s_cells[mw,2,i] < 0: + q[:,iw+1] -= (arrival_times[i+1] - arrival_times[i]) / ((1 - alpha) * dx) * (wave_cells[:,mw,2,i]) + if s_cells[mw,1,i] < 0: + q[:,iw-1] -= (dt - arrival_times[i]) / dx * ( max(0, -s_cells[mw,1,i] * (dt - arrival_times[i]) - alpha * dx) / (-s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) + q[:,iw] -= (dt - arrival_times[i]) / (alpha * dx) * ( min(-s_cells[mw,1, i] * (dt - arrival_times[i]), alpha * dx) / (-s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) + if s_cells[mw,1,i] > 0: + q[:,iw+1] -= (dt - arrival_times[i]) / ((1 - alpha) * dx) * ( min(s_cells[mw,1, i] * (dt - arrival_times[i]), (1 - alpha) * dx) / (s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) + q[:,iw+2] -= (dt - arrival_times[i]) / dx * ( max(0, s_cells[mw,1,i] * (dt - arrival_times[i]) - (1- alpha) * dx) / (s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) + + + # Compute maximum wave speed + # add additional conditions for h-box + cfl = 0.0 + if 'method' not in state.problem_data: + for mw in range(wave.shape[1]): + smax1 = np.max(dtdx[LL:UL]*s[mw,LL-1:UL-1]) + smax2 = np.max(-dtdx[LL-1:UL-1]*s[mw,LL-1:UL-1]) + cfl = max(cfl,smax1,smax2) + elif state.problem_data['method'] == 'h_box': + for mw in range(wave.shape[1]): + smax1 = np.max(dtdx_hbox[LL:UL]*s[mw,LL-1:UL-1]) + smax2 = np.max(-dtdx_hbox[LL-1:UL-1]*s[mw,LL-1:UL-1]) + cfl = max(cfl,smax1,smax2) + elif state.problem_data['method'] == 'h_box_wave': + for mw in range(wave.shape[1]): + smax1 = np.max(dtdx_hbox[LL:UL]*s[mw,LL-1:UL-1]) + smax2 = np.max(-dtdx_hbox[LL-1:UL-1]*s[mw,LL-1:UL-1]) + cfl = max(cfl,smax1,smax2) + + + # If we are doing slope limiting we have more work to do + if self.order == 2: + # Initialize flux corrections + f = np.zeros( (num_eqn,grid.num_cells[0] + 2*self.num_ghost) ) + + # Apply Limiters to waves + if (limiter > 0).any(): + wave = tvd.limit(state.num_eqn,wave,s,limiter,dtdx) + + # Compute correction fluxes for second order q_{xx} terms + dtdxave = 0.5 * (dtdx[LL-1:UL-1] + dtdx[LL:UL]) + if self.fwave: + for mw in range(wave.shape[1]): + sabs = np.abs(s[mw,LL-1:UL-1]) + om = 1.0 - sabs*dtdxave[:UL-LL] + ssign = np.sign(s[mw,LL-1:UL-1]) + for m in range(num_eqn): + f[m,LL:UL] += 0.5 * ssign * om * wave[m,mw,LL-1:UL-1] + else: + for mw in range(wave.shape[1]): + sabs = np.abs(s[mw,LL-1:UL-1]) + om = 1.0 - sabs*dtdxave[:UL-LL] + for m in range(num_eqn): + f[m,LL:UL] += 0.5 * sabs * om * wave[m,mw,LL-1:UL-1] + + # Update q by differencing correction fluxes + for m in range(num_eqn): + q[m,LL:UL-1] -= dtdx[LL:UL-1] * (f[m,LL+1:UL] - f[m,LL:UL-1]) + + + else: raise Exception("Unrecognized kernel_language; choose 'Fortran' or 'Python'") + + self.cfl.update_global_max(cfl) + state.set_q_from_qbc(num_ghost,self.qbc) + if state.num_aux > 0: + state.set_aux_from_auxbc(num_ghost,self.auxbc) + + +# ============================================================================ +# ClawPack 2d Solver Class +# ============================================================================ +class ClawSolver2D(ClawSolver): + r""" + 2D Classic (Clawpack) solver. + + Solve using the wave propagation algorithms of Randy LeVeque's + Clawpack code (www.clawpack.org). + + In addition to the attributes of ClawSolver1D, ClawSolver2D + also has the following options: + + .. attribute:: dimensional_split + + If True, use dimensional splitting (Godunov splitting). + Dimensional splitting with Strang splitting is not supported + at present but could easily be enabled if necessary. + If False, use unsplit Clawpack algorithms, possibly including + transverse Riemann solves. + + .. attribute:: transverse_waves + + If dimensional_split is True, this option has no effect. If + dimensional_split is False, then transverse_waves should be one of + the following values: + + ClawSolver2D.no_trans: Transverse Riemann solver + not used. The stable CFL for this algorithm is 0.5. Not recommended. + + ClawSolver2D.trans_inc: Transverse increment waves are computed + and propagated. + + ClawSolver2D.trans_cor: Transverse increment waves and transverse + correction waves are computed and propagated. + + Note that only the fortran routines are supported for now in 2D. + """ + + __doc__ += add_parent_doc(ClawSolver) + + no_trans = 0 + trans_inc = 1 + trans_cor = 2 + + def __init__(self,riemann_solver=None, claw_package=None): + r""" + Create 2d Clawpack solver + + See :class:`ClawSolver2D` for more info. + """ + self.dimensional_split = True + self.transverse_waves = self.trans_inc + + self.num_dim = 2 + self.reflect_index = [1,2] + + self.aux1 = None + self.aux2 = None + self.aux3 = None + self.work = None + + super(ClawSolver2D,self).__init__(riemann_solver, claw_package) + + def _check_cfl_settings(self): + if (not self.dimensional_split) and (self.transverse_waves==0): + cfl_recommended = 0.5 + else: + cfl_recommended = 1.0 + + if self.cfl_max > cfl_recommended: + import warnings + warnings.warn('cfl_max is set higher than the recommended value of %s' % cfl_recommended) + warnings.warn(str(self.cfl_desired)) + + + def _allocate_workspace(self,solution): + r""" + Pack parameters into format recognized by Clawpack (Fortran) code. + + Sets the method array and the cparam common block for the Riemann solver. + """ + import numpy as np + + state = solution.state + + num_eqn,num_aux,num_waves,num_ghost,aux = state.num_eqn,state.num_aux,self.num_waves,self.num_ghost,state.aux + + #The following is a hack to work around an issue + #with f2py. It involves wastefully allocating three arrays. + #f2py seems not able to handle multiple zero-size arrays being passed. + # it appears the bug is related to f2py/src/fortranobject.c line 841. + if aux is None: num_aux=1 + + grid = state.grid + maxmx,maxmy = grid.num_cells[0],grid.num_cells[1] + maxm = max(maxmx, maxmy) + + # These work arrays really ought to live inside a fortran module + # as is done for sharpclaw + self.aux1 = np.empty((num_aux,maxm+2*num_ghost),order='F') + self.aux2 = np.empty((num_aux,maxm+2*num_ghost),order='F') + self.aux3 = np.empty((num_aux,maxm+2*num_ghost),order='F') + mwork = (maxm+2*num_ghost) * (5*num_eqn + num_waves + num_eqn*num_waves) + self.work = np.empty((mwork),order='F') + + + # ========== Hyperbolic Step ===================================== + def step_hyperbolic(self,solution): + r""" + Take a step on the homogeneous hyperbolic system using the Clawpack + algorithm. + + Clawpack is based on the Lax-Wendroff method, combined with Riemann + solvers and TVD limiters applied to waves. + """ + if(self.kernel_language == 'Fortran'): + state = solution.states[0] + grid = state.grid + dx,dy = grid.delta + mx,my = grid.num_cells + maxm = max(mx,my) + + self._apply_bcs(state) + qold = self.qbc.copy('F') + + rpn2 = self.rp.rpn2._cpointer + + if (self.dimensional_split) or (self.transverse_waves==0): + rpt2 = rpn2 # dummy value; it won't be called + else: + rpt2 = self.rp.rpt2._cpointer + + if self.dimensional_split: + #Right now only Godunov-dimensional-splitting is implemented. + #Strang-dimensional-splitting could be added following dimsp2.f in Clawpack. + + self.qbc, cfl_x = self.fmod.step2ds(maxm,self.num_ghost,mx,my, \ + qold,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,1,self.fwave,rpn2,rpt2) + + self.qbc, cfl_y = self.fmod.step2ds(maxm,self.num_ghost,mx,my, \ + self.qbc,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,2,self.fwave,rpn2,rpt2) + + cfl = max(cfl_x,cfl_y) + + else: + + self.qbc, cfl = self.fmod.step2(maxm,self.num_ghost,mx,my, \ + qold,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,self.fwave,rpn2,rpt2) + + self.cfl.update_global_max(cfl) + state.set_q_from_qbc(self.num_ghost,self.qbc) + if state.num_aux > 0: + state.set_aux_from_auxbc(self.num_ghost,self.auxbc) + + else: + raise NotImplementedError("No python implementation for step_hyperbolic in 2D.") + +# ============================================================================ +# ClawPack 3d Solver Class +# ============================================================================ +class ClawSolver3D(ClawSolver): + r""" + 3D Classic (Clawpack) solver. + + Solve using the wave propagation algorithms of Randy LeVeque's + Clawpack code (www.clawpack.org). + + In addition to the attributes of ClawSolver, ClawSolver3D + also has the following options: + + .. attribute:: dimensional_split + + If True, use dimensional splitting (Godunov splitting). + Dimensional splitting with Strang splitting is not supported + at present but could easily be enabled if necessary. + If False, use unsplit Clawpack algorithms, possibly including + transverse Riemann solves. + + .. attribute:: transverse_waves + + If dimensional_split is True, this option has no effect. If + dim_plit is False, then transverse_waves should be one of + the following values: + + ClawSolver3D.no_trans: Transverse Riemann solver + not used. The stable CFL for this algorithm is 0.5. Not recommended. + + ClawSolver3D.trans_inc: Transverse increment waves are computed + and propagated. + + ClawSolver3D.trans_cor: Transverse increment waves and transverse + correction waves are computed and propagated. + + Note that only Fortran routines are supported for now in 3D -- + there is no pure-python version. + """ + + __doc__ += add_parent_doc(ClawSolver) + + no_trans = 0 + trans_inc = 11 + trans_cor = 22 + + def __init__(self, riemann_solver=None, claw_package=None): + r""" + Create 3d Clawpack solver + + See :class:`ClawSolver3D` for more info. + """ + # Add the functions as required attributes + self.dimensional_split = True + self.transverse_waves = self.trans_cor + + self.num_dim = 3 + self.reflect_index = [1,2,3] + + self.aux1 = None + self.aux2 = None + self.aux3 = None + self.work = None + + super(ClawSolver3D,self).__init__(riemann_solver, claw_package) + + # ========== Setup routine ============================= + def _allocate_workspace(self,solution): + r""" + Allocate auxN and work arrays for use in Fortran subroutines. + """ + import numpy as np + + state = solution.states[0] + + num_eqn,num_aux,num_waves,num_ghost,aux = state.num_eqn,state.num_aux,self.num_waves,self.num_ghost,state.aux + + #The following is a hack to work around an issue + #with f2py. It involves wastefully allocating three arrays. + #f2py seems not able to handle multiple zero-size arrays being passed. + # it appears the bug is related to f2py/src/fortranobject.c line 841. + if(aux is None): num_aux=1 + + grid = state.grid + maxmx,maxmy,maxmz = grid.num_cells[0],grid.num_cells[1],grid.num_cells[2] + maxm = max(maxmx, maxmy, maxmz) + + # These work arrays really ought to live inside a fortran module + # as is done for sharpclaw + self.aux1 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') + self.aux2 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') + self.aux3 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') + mwork = (maxm+2*num_ghost) * (31*num_eqn + num_waves + num_eqn*num_waves) + self.work = np.empty((mwork),order='F') + + + # ========== Hyperbolic Step ===================================== + def step_hyperbolic(self,solution): + r""" + Take a step on the homogeneous hyperbolic system using the Clawpack + algorithm. + + Clawpack is based on the Lax-Wendroff method, combined with Riemann + solvers and TVD limiters applied to waves. + """ + if(self.kernel_language == 'Fortran'): + state = solution.states[0] + grid = state.grid + dx,dy,dz = grid.delta + mx,my,mz = grid.num_cells + maxm = max(mx,my,mz) + + self._apply_bcs(state) + qnew = self.qbc + qold = qnew.copy('F') + + rpn3 = self.rp.rpn3._cpointer + + if (self.dimensional_split) or (self.transverse_waves==0): + rpt3 = rpn3 # dummy value; it won't be called + rptt3 = rpn3 # dummy value; it won't be called + else: + rpt3 = self.rp.rpt3._cpointer + rptt3 = self.rp.rptt3._cpointer + + if self.dimensional_split: + #Right now only Godunov-dimensional-splitting is implemented. + #Strang-dimensional-splitting could be added following dimsp3.f in Clawpack. + + q, cfl_x = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ + qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,1,self.fwave,rpn3,rpt3,rptt3) + + q, cfl_y = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ + q,q,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,2,self.fwave,rpn3,rpt3,rptt3) + + q, cfl_z = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ + q,q,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,3,self.fwave,rpn3,rpt3,rptt3) + + cfl = max(cfl_x,cfl_y,cfl_z) + + else: + + q, cfl = self.fmod.step3(maxm,self.num_ghost,mx,my,mz, \ + qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ + self.aux1,self.aux2,self.aux3,self.work,self.fwave,rpn3,rpt3,rptt3) + + self.cfl.update_global_max(cfl) + state.set_q_from_qbc(self.num_ghost,self.qbc) + if state.num_aux > 0: + state.set_aux_from_auxbc(self.num_ghost,self.auxbc) + + else: + raise NotImplementedError("No python implementation for step_hyperbolic in 3D.") diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 02b90f5bd..d1b5c5a7f 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -31,7 +31,7 @@ class SharpClawSolver(Solver): r""" Superclass for all SharpClawND solvers. - Implements Runge-Kutta time stepping and the basic form of a + Implements Runge-Kutta time stepping and the basic form of a semi-discrete step (the dq() function). If another method-of-lines solver is implemented in the future, it should be based on this class, which then ought to be renamed to something like "MOLSolver". @@ -93,11 +93,11 @@ class SharpClawSolver(Solver): Whether the auxiliary array is time dependent. ``Default = False`` - + .. attribute:: kernel_language Specifies whether to use wrapped Fortran routines ('Fortran') - or pure Python ('Python'). + or pure Python ('Python'). ``Default = 'Fortran'``. .. attribute:: num_ghost @@ -106,9 +106,9 @@ class SharpClawSolver(Solver): ``Default = 3`` .. attribute:: fwave - - Whether to split the flux jump (rather than the jump in Q) into waves; - requires that the Riemann solver performs the splitting. + + Whether to split the flux jump (rather than the jump in Q) into waves; + requires that the Riemann solver performs the splitting. ``Default = False`` .. attribute:: cfl_desired @@ -123,7 +123,7 @@ class SharpClawSolver(Solver): .. attribute:: dq_src - Whether a source term is present. If it is present the function that + Whether a source term is present. If it is present the function that computes its contribution must be provided. ``Default = None`` @@ -178,7 +178,7 @@ def __init__(self,riemann_solver=None,claw_package=None): self._method = None self._registers = None self.dq_dt = None - self.dt_old = None + self.dt_old = None # Used only if time integrator is 'RK' self.a = None @@ -200,7 +200,7 @@ def __init__(self,riemann_solver=None,claw_package=None): # Call general initialization function super(SharpClawSolver,self).__init__(riemann_solver,claw_package) - + def setup(self,solution): """ Allocate RK stage arrays or previous step solutions and fortran routine work arrays. @@ -226,19 +226,19 @@ def setup(self,solution): assert 3 <= self.lmm_steps, \ 'Must set solver.lmm_steps greater than 2 for 2nd order SSPLMM integrator.' self.sspcoeff = (self.lmm_steps - 2.)/(self.lmm_steps - 1.) - # The choice of cfl_desired and cfl_max is intended for LMM with many steps (up to 20). + # The choice of cfl_desired and cfl_max is intended for LMM with many steps (up to 20). # If more steps are chosen the solution may not be accurate enough. # For a smaller number of steps, higher values of cfl_desired and cfl_max can be used. if self.cfl_max is None: - self.cfl_desired = 0.14*self.sspcoeff - self.cfl_max = 0.15*self.sspcoeff + self.cfl_desired = 0.14*self.sspcoeff + self.cfl_max = 0.15*self.sspcoeff elif self.time_integrator == 'SSPLMMk3': assert 4 <= self.lmm_steps <= 5, \ 'Must set solver.lmm_steps equal to 4 or 5 for 3rd order SSPLMM integrator.' self.sspcoeff = (self.lmm_steps - 3.)/(self.lmm_steps - 1.) if self.cfl_max is None: - self.cfl_desired = 0.48*self.sspcoeff - self.cfl_max = 0.5*self.sspcoeff + self.cfl_desired = 0.48*self.sspcoeff + self.cfl_max = 0.5*self.sspcoeff else: if self.cfl_max is None: self.cfl_desired = self._cfl_default[self.time_integrator][0] @@ -252,7 +252,7 @@ def setup(self,solution): self._set_mthlim() state = solution.states[0] - + if self.kernel_language=='Fortran': if self.fmod is None: so_name = 'clawpack.pyclaw.sharpclaw.sharpclaw'+str(self.num_dim) @@ -431,7 +431,7 @@ def ssp104(self,state): def update_saved_values(self,state,step_index): - r""" + r""" Updates lists of saved function evaluations, solution values, dt and dtFE for LMMs. For 3rd-order SSPLMM additional conditions are checked if self.check_lmm_cond is set to True. If these conditions are violated, the step is rejected. @@ -496,7 +496,7 @@ def check_3rd_ord_cond(self,state,step_index,dtFE): r""" This routine checks the additional conditions for the 3rd-order SSPLMMs. This is a posteriori check after a step is accepted. - In particular, there is a condition on the step size for the starting values and + In particular, there is a condition on the step size for the starting values and a condition on the ratio of forward Euler step sizes at very step. If the conditions are violated we muct retrieve the previous solution and discard that step; otherwise the step is accepted. @@ -523,7 +523,7 @@ def _set_mthlim(self): if len(self._mthlim)!=self.num_waves: raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves') - + def dq(self,state): """ Evaluate dq/dt * (delta t) @@ -540,7 +540,7 @@ def dq(self,state): def dq_hyperbolic(self,state): raise NotImplementedError('You must subclass SharpClawSolver.') - + def dqdt(self,state): """ Evaluate dq/dt. This routine is used for implicit time stepping. @@ -589,12 +589,12 @@ def _allocate_registers(self,solution): This routine is only used by method-of-lines solvers (SharpClaw), not by the Classic solvers. It allocates additional State objects - to store the intermediate stages used by Runge--Kutta and Multistep + to store the intermediate stages used by Runge--Kutta and Multistep time integrators. If we create a MethodOfLinesSolver subclass, this should be moved there. """ - # Generally the number of registers for the starting method should be at most + # Generally the number of registers for the starting method should be at most # equal to the number of registers of the LMM if self.time_integrator == 'Euler': nregisters=0 elif self.time_integrator == 'SSP33': nregisters=1 @@ -605,7 +605,7 @@ def _allocate_registers(self,solution): elif self.time_integrator == 'LMM': nregisters=len(self.alpha) else: raise Exception('Unrecognized time integrator: '+self.time_integrator) - + state = solution.states[0] # Use the same class constructor as the solution for the Runge Kutta stages self._registers = [] @@ -619,7 +619,7 @@ def accept_reject_step(self,state): Decide whether to accept or not the current step. For Runge-Kutta methods the step is accepted if cfl <= cfl_max. For SSPLMM32 the choice of step-size guarantees the cfl condition is satisfied for the steps the LMM - is used. Hence, we need to check the cfl condition only for the starting steps. + is used. Hence, we need to check the cfl condition only for the starting steps. """ accept_step = True cfl = self.cfl.get_cached_max() @@ -650,7 +650,7 @@ def get_dt_new(self): Set size of next step depending on the time integrator and whether or not the current step was accepted. """ - self.dt_old = self.dt + self.dt_old = self.dt cfl = self.cfl.get_cached_max() if 'SSPLMM' in self.time_integrator: @@ -686,17 +686,17 @@ class SharpClawSolver1D(SharpClawSolver): # ======================================================================== """ SharpClaw solver for one-dimensional problems. - + Used to solve 1D hyperbolic systems using the SharpClaw algorithms, which are based on WENO reconstruction and Runge-Kutta time stepping. """ __doc__ += add_parent_doc(SharpClawSolver) - + def __init__(self,riemann_solver=None,claw_package=None): r""" See :class:`SharpClawSolver1D` for more info. - """ + """ self.num_dim = 1 self.reflect_index = [1] super(SharpClawSolver1D,self).__init__(riemann_solver,claw_package) @@ -714,11 +714,11 @@ def dq_hyperbolic(self,state): 0 1 2 3 4 mx+num_ghost-2 mx+num_ghost mx+num_ghost+2 | mx+num_ghost-1 | mx+num_ghost+1 | | | | | ... | | | | | - 0 1 | 2 3 mx+num_ghost-2 |mx+num_ghost + 0 1 | 2 3 mx+num_ghost-2 |mx+num_ghost mx+num_ghost-1 mx+num_ghost+1 The top indices represent the values that are located on the grid - cell boundaries such as waves, s and other Riemann problem values, + cell boundaries such as waves, s and other Riemann problem values, the bottom for the cell centered values such as q. In particular the ith grid cell boundary has the following related information:: @@ -731,11 +731,11 @@ def dq_hyperbolic(self,state): values are in the cell. """ - + import numpy as np self._apply_bcs(state) - q = self.qbc + q = self.qbc grid = state.grid mx = grid.num_cells[0] @@ -752,17 +752,16 @@ def dq_hyperbolic(self,state): dq,cfl=self.fmod.flux1(q,self.auxbc,self.dt,state.t,ixy,mx,self.num_ghost,mx,rp1,tfluct1) elif self.kernel_language=='Python': - + aux=self.auxbc dtdx = np.zeros( (mx+2*self.num_ghost) ,order='F') dq = np.zeros( (state.num_eqn,mx+2*self.num_ghost) ,order='F') # Find local value for dt/dx if state.index_capa>=0: - dtdx = self.dt / (grid.delta[0] * state.aux[state.index_capa,:]) + dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) else: dtdx += self.dt/grid.delta[0] - - aux=self.auxbc + if aux.shape[0]>0: aux_l=aux[:,:-1] aux_r=aux[:,1: ] @@ -813,12 +812,12 @@ def dq_hyperbolic(self,state): dq[m,LL:UL] = -dtdx[LL:UL]*(amdq[m,LL:UL] + apdq[m,LL-1:UL-1] \ + apdq2[m,LL:UL] + amdq2[m,LL:UL]) - else: + else: raise Exception('Unrecognized value of solver.kernel_language.') self.cfl.update_global_max(cfl) return dq[:,self.num_ghost:-self.num_ghost] - + # ======================================================================== class SharpClawSolver2D(SharpClawSolver): @@ -831,9 +830,9 @@ class SharpClawSolver2D(SharpClawSolver): def __init__(self,riemann_solver=None,claw_package=None): r""" Create 2D SharpClaw solver - + See :class:`SharpClawSolver2D` for more info. - """ + """ self.num_dim = 2 self.reflect_index = [1,2] @@ -851,11 +850,11 @@ def dq_hyperbolic(self,state): 0 1 2 3 4 mx+num_ghost-2 mx+num_ghost mx+num_ghost+2 | mx+num_ghost-1 | mx+num_ghost+1 | | | | | ... | | | | | - 0 1 | 2 3 mx+num_ghost-2 |mx+num_ghost + 0 1 | 2 3 mx+num_ghost-2 |mx+num_ghost mx+num_ghost-1 mx+num_ghost+1 The top indices represent the values that are located on the grid - cell boundaries such as waves, s and other Riemann problem values, + cell boundaries such as waves, s and other Riemann problem values, the bottom for the cell centered values such as q. In particular the ith grid cell boundary has the following related information:: @@ -869,7 +868,7 @@ def dq_hyperbolic(self,state): """ self._apply_bcs(state) - q = self.qbc + q = self.qbc grid = state.grid @@ -903,9 +902,9 @@ class SharpClawSolver3D(SharpClawSolver): def __init__(self,riemann_solver=None,claw_package=None): r""" Create 3D SharpClaw solver - + See :class:`SharpClawSolver3D` for more info. - """ + """ self.num_dim = 3 self.reflect_index = [1,2,3] @@ -935,11 +934,11 @@ def dq_hyperbolic(self,state): 0 1 2 3 4 mx+num_ghost-2 mx+num_ghost mx+num_ghost+2 | mx+num_ghost-1 | mx+num_ghost+1 | | | | | ... | | | | | - 0 1 | 2 3 mx+num_ghost-2 |mx+num_ghost + 0 1 | 2 3 mx+num_ghost-2 |mx+num_ghost mx+num_ghost-1 mx+num_ghost+1 The top indices represent the values that are located on the grid - cell boundaries such as waves, s and other Riemann problem values, + cell boundaries such as waves, s and other Riemann problem values, the bottom for the cell centered values such as q. In particular the ith grid cell boundary has the following related information:: @@ -953,7 +952,7 @@ def dq_hyperbolic(self,state): """ self._apply_bcs(state) - q = self.qbc + q = self.qbc grid = state.grid @@ -965,7 +964,7 @@ def dq_hyperbolic(self,state): if self.kernel_language=='Fortran': rpn3 = self.rp.rpn3._cpointer - if self.tfluct_solver: + if self.tfluct_solver: tfluct3 = self.tfluct.tfluct3._cpointer else: tfluct3 = self.tfluct From ec122e8dbf47ec614f9b402e6cc238ae4e5f0cf0 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Fri, 18 Oct 2019 22:53:25 -0400 Subject: [PATCH 41/47] Delete solverjli.py --- src/pyclaw/classic/solverjli.py | 866 -------------------------------- 1 file changed, 866 deletions(-) delete mode 100644 src/pyclaw/classic/solverjli.py diff --git a/src/pyclaw/classic/solverjli.py b/src/pyclaw/classic/solverjli.py deleted file mode 100644 index 7dfd80c57..000000000 --- a/src/pyclaw/classic/solverjli.py +++ /dev/null @@ -1,866 +0,0 @@ -r""" -Module containing the classic Clawpack solvers. - -This module contains the pure and wrapped classic clawpack solvers. All -clawpack solvers inherit from the :class:`ClawSolver` superclass which in turn -inherits from the :class:`~pyclaw.solver.Solver` superclass. These -are both pure virtual classes; the only solver classes that should be instantiated -are the dimension-specific ones, :class:`ClawSolver1D` and :class:`ClawSolver2D`. -""" - -from clawpack.pyclaw.util import add_parent_doc -from clawpack.pyclaw.solver import Solver -from clawpack.pyclaw.limiters import tvd - -# ============================================================================ -# Generic Clawpack solver class -# ============================================================================ -class ClawSolver(Solver): - r""" - Generic classic Clawpack solver - - All Clawpack solvers inherit from this base class. - - .. attribute:: mthlim - - Limiter(s) to be used. Specified either as one value or a list. - If one value, the specified limiter is used for all wave families. - If a list, the specified values indicate which limiter to apply to - each wave family. Take a look at pyclaw.limiters.tvd for an enumeration. - ``Default = limiters.tvd.minmod`` - - .. attribute:: order - - Order of the solver, either 1 for first order (i.e., Godunov's method) - or 2 for second order (Lax-Wendroff-LeVeque). - ``Default = 2`` - - .. attribute:: source_split - - Which source splitting method to use: 1 for first - order Godunov splitting and 2 for second order Strang splitting. - ``Default = 1`` - - .. attribute:: fwave - - Whether to split the flux jump (rather than the jump in Q) into waves; - requires that the Riemann solver performs the splitting. - ``Default = False`` - - .. attribute:: step_source - - Handle for function that evaluates the source term. - The required signature for this function is: - - def step_source(solver,state,dt) - - .. attribute:: kernel_language - - Specifies whether to use wrapped Fortran routines ('Fortran') - or pure Python ('Python'). ``Default = 'Fortran'``. - - .. attribute:: verbosity - - The level of detail of logged messages from the Fortran solver. - ``Default = 0``. - - """ - - # ========== Generic Init Routine ======================================== - def __init__(self,riemann_solver=None,claw_package=None): - r""" - See :class:`ClawSolver` for full documentation. - - Output: - - (:class:`ClawSolver`) - Initialized clawpack solver - """ - self.num_ghost = 2 - self.limiters = tvd.vanleer - self.order = 2 - self.source_split = 1 - self.fwave = False - self.step_source = None - self.kernel_language = 'Fortran' - self.verbosity = 0 - self.cfl_max = 1.0 - self.cfl_desired = 0.9 - self._mthlim = self.limiters - self._method = None - self.dt_old = None - - # Call general initialization function - super(ClawSolver,self).__init__(riemann_solver,claw_package) - - # ========== Time stepping routines ====================================== - def step(self,solution,take_one_step,tstart,tend): - r""" - Evolve solution one time step - - The elements of the algorithm for taking one step are: - - 1. Pick a step size as specified by the base solver attribute :func:`get_dt` - - 2. A half step on the source term :func:`step_source` if Strang splitting is - being used (:attr:`source_split` = 2) - - 3. A step on the homogeneous problem :math:`q_t + f(q)_x = 0` is taken - - 4. A second half step or a full step is taken on the source term - :func:`step_source` depending on whether Strang splitting was used - (:attr:`source_split` = 2) or Godunov splitting - (:attr:`source_split` = 1) - - This routine is called from the method evolve_to_time defined in the - pyclaw.solver.Solver superclass. - - :Input: - - *solution* - (:class:`~pyclaw.solution.Solution`) solution to be evolved - - :Output: - - (bool) - True if full step succeeded, False otherwise - """ - self.get_dt(solution.t,tstart,tend,take_one_step) - self.cfl.set_global_max(0.) - - if self.source_split == 2 and self.step_source is not None: - self.step_source(self,solution.states[0],self.dt/2.0) - - self.step_hyperbolic(solution) - - # Check here if the CFL condition is satisfied. - # If not, return # immediately to evolve_to_time and let it deal with - # picking a new step size (dt). - if self.cfl.get_cached_max() >= self.cfl_max: - return False - - if self.step_source is not None: - # Strang splitting - if self.source_split == 2: - self.step_source(self,solution.states[0],self.dt/2.0) - - # Godunov Splitting - if self.source_split == 1: - self.step_source(self,solution.states[0],self.dt) - - return True - - def _check_cfl_settings(self): - pass - - def _allocate_workspace(self,solution): - pass - - def step_hyperbolic(self,solution): - r""" - Take one homogeneous step on the solution. - - This is a dummy routine and must be overridden. - """ - raise Exception("Dummy routine, please override!") - - def _set_mthlim(self): - r""" - Convenience routine to convert users limiter specification to - the format understood by the Fortran code (i.e., a list of length num_waves). - """ - #self._mthlim = self.limiters - #if not isinstance(self.limiters,list): self._mthlim=[self._mthlim] - # if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves - # if len(self._mthlim)!=self.num_waves: - # raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves') - - def _set_method(self,state): - r""" - Set values of the solver._method array required by the Fortran code. - These are algorithmic parameters. - """ - import numpy as np - #We ought to put method and many other things in a Fortran - #module and set the fortran variables directly here. - self._method =np.empty(7, dtype=int,order='F') - self._method[0] = self.dt_variable - self._method[1] = self.order - if self.num_dim==1: - self._method[2] = 0 # Not used in 1D - elif self.dimensional_split: - self._method[2] = -1 # First-order dimensional splitting - else: - self._method[2] = self.transverse_waves - self._method[3] = self.verbosity - self._method[4] = 0 # Not used for PyClaw (would be self.source_split) - self._method[5] = state.index_capa + 1 - self._method[6] = state.num_aux - - def setup(self,solution): - r""" - Perform essential solver setup. This routine must be called before - solver.step() may be called. - """ - # This is a hack to deal with the fact that petsc4py - # doesn't allow us to change the stencil_width (num_ghost) - solution.state.set_num_ghost(self.num_ghost) - # End hack - - self._check_cfl_settings() - - self._set_mthlim() - if(self.kernel_language == 'Fortran'): - if self.fmod is None: - so_name = 'clawpack.pyclaw.classic.classic'+str(self.num_dim) - self.fmod = __import__(so_name,fromlist=['clawpack.pyclaw.classic']) - self._set_fortran_parameters(solution) - self._allocate_workspace(solution) - elif self.num_dim>1: - raise Exception('Only Fortran kernels are supported in multi-D.') - - self._allocate_bc_arrays(solution.states[0]) - - super(ClawSolver,self).setup(solution) - - - def _set_fortran_parameters(self,solution): - r""" - Pack parameters into format recognized by Clawpack (Fortran) code. - - Sets the solver._method array and the cparam common block for the Riemann solver. - """ - self._set_method(solution.state) - # The reload here is necessary because otherwise the common block - # cparam in the Riemann solver doesn't get flushed between running - # different tests in a single Python session. - reload(self.fmod) - solution.state.set_cparam(self.fmod) - solution.state.set_cparam(self.rp) - - def __del__(self): - r""" - Delete Fortran objects, which otherwise tend to persist in Python sessions. - """ - if(self.kernel_language == 'Fortran'): - del self.fmod - - super(ClawSolver,self).__del__() - - -# ============================================================================ -# ClawPack 1d Solver Class -# ============================================================================ -class ClawSolver1D(ClawSolver): - r""" - Clawpack evolution routine in 1D - - This class represents the 1d clawpack solver on a single grid. Note that - there are routines here for interfacing with the fortran time stepping - routines and the Python time stepping routines. The ones used are - dependent on the argument given to the initialization of the solver - (defaults to python). - - """ - - __doc__ += add_parent_doc(ClawSolver) - - def __init__(self, riemann_solver=None, claw_package=None): - r""" - Create 1d Clawpack solver - - Output: - - (:class:`ClawSolver1D`) - Initialized 1d clawpack solver - - See :class:`ClawSolver1D` for more info. - """ - self.num_dim = 1 - self.reflect_index = [1] - - super(ClawSolver1D,self).__init__(riemann_solver, claw_package) - - - # ========== Homogeneous Step ===================================== - def step_hyperbolic(self,solution): - r""" - Take one time step on the homogeneous hyperbolic system. - - :Input: - - *solution* - (:class:`~pyclaw.solution.Solution`) Solution that - will be evolved - """ - import numpy as np - - state = solution.states[0] - grid = state.grid - - - self._apply_bcs(state) - - num_eqn,num_ghost = state.num_eqn,self.num_ghost - - if(self.kernel_language == 'Fortran'): - mx = grid.num_cells[0] - dx,dt = grid.delta[0],self.dt - dtdx = np.zeros( (mx+2*num_ghost) ) + dt/dx - rp1 = self.rp.rp1._cpointer - - self.qbc,cfl = self.fmod.step1(num_ghost,mx,self.qbc,self.auxbc,dx,dt,self._method,self._mthlim,self.fwave,rp1) - - elif(self.kernel_language == 'Python'): - q = self.qbc - aux = self.auxbc - # Limiter to use in the pth family - limiter = np.array(self._mthlim,ndmin=1) - - dtdx = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - - # Find local value for dt/dx - if 'method' not in state.problem_data: - if state.index_capa>=0: - dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) - else: - dtdx += self.dt / grid.delta[0] - elif state.problem_data['method'] == 'h_box': - xpxc = state.problem_data['xpxc'] - dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) - dtdx_hbox = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - dtdx_hbox += self.dt / (grid.delta[0] * xpxc) - elif state.problem_data['method'] == 'h_box_wave': - xpxc = state.problem_data['xpxc'] - dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) - dtdx_hbox = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - dtdx_hbox += self.dt / (grid.delta[0] * xpxc) - - - - # Solve Riemann problem at each interface - q_l=q[:,:-1].copy() - q_r=q[:,1:].copy() - if state.aux is not None: - aux_l=aux[:,:-1].copy() - aux_r=aux[:,1:].copy() - else: - aux_l = None - aux_r = None - - if 'method' not in state.problem_data: - # normal case - wave,s,amdq,apdq = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - - - # Update loop limits, these are the limits for the Riemann solver - # locations, which then update a grid cell value - # We include the Riemann problem just outside of the grid so we can - # do proper limiting at the grid edges - # LL | | UL - # | LL | | | | ... | | | UL | | - # | | - - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - - # Update q for Godunov update - for m in range(num_eqn): - q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] - q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] - - elif state.problem_data['method'] == 'h_box': - # # add corrections - wave,s,amdq,apdq,f_corr_l,f_corr_r = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - - # Update q for Godunov update - for m in range(num_eqn): - q[m,LL:UL] -= dtdx[LL:UL]*(apdq[m,LL-1:UL-1] - f_corr_r[m,LL-1:UL-1]) - q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*(amdq[m,LL-1:UL-1] + f_corr_l[m,LL-1:UL-1]) - - elif state.problem_data['method'] == 'h_box_wave': - # # add corrections - state.problem_data['arrival_state'] = False - - - wave,s,amdq,apdq,q_hbox_initial,aux_hbox = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - - # Update q for Godunov update - iw = state.problem_data['wall_position'] + self.num_ghost - 1 - q_last = q[:,iw:iw+2].copy() - - for m in range(num_eqn): - q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] - q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] - - # check the arrivals - q[:,iw:iw+2] = q_last[:,:] # reset the wall cells - dt = self.dt - num_waves = self.num_waves - dx = grid.delta[0] * xpxc - alpha = state.problem_data['fraction'] - arrival_times = np.array([0.0]) - for mw in range(num_waves): - if (s[mw,iw-1] > 0 and (s[mw,iw-1] * dt > alpha * dx)): - arrival_times = np.append(arrival_times, alpha * dx / s[mw,iw-1]) - if (s[mw,iw+1] < 0 and ( (-s[mw,iw+1]) * dt > (1 - alpha) * dx ) ): - arrival_times = np.append(arrival_times, -(1 - alpha) * dx / s[mw,iw+1]) - arrival_times.sort() - - n_arrival_times = len(arrival_times) - if n_arrival_times == 1 : - state.problem_data['arrival_state'] = False - else: - state.problem_data['arrival_state'] = True - s_cells = np.zeros((num_waves, 3, n_arrival_times)) - s_cells[:,:,0] = s[:, iw-1:iw+2].copy() - wave_cells = np.zeros((num_eqn, num_waves, 3, n_arrival_times)) - wave_cells[:,:,:,0] = wave[:,:,iw-1:iw+2].copy() - - - if state.problem_data['arrival_state'] == False: - q[:,iw] -= dt/(alpha * dx) * apdq[:,iw-1] - q[:,iw+1] -= dt/((1 - alpha)*dx) * amdq[:,iw+1] - for mw in range(num_waves): - if (s[mw,iw] < 0): - q[:,iw-1] -= dt/dx * ( max(0, -s[mw,iw] * dt - alpha * dx) / (-s[mw,iw] * dt) * wave[:,mw,iw] ) - q[:,iw] -= dt/(alpha * dx) * (min(-s[mw,iw] * dt, alpha * dx) / (-s[mw,iw] * dt) * wave[:,mw,iw] ) - elif (s[mw,iw] > 0): - q[:,iw+1] -= dt/((1 - alpha)*dx) * (min(s[mw,iw] * dt, (1 - alpha) * dx) / (s[mw,iw] * dt) * wave[:,mw,iw] ) - q[:,iw+2] -= dt/dx * ( (max(0, s[mw,iw] * dt - (1 - alpha) * dx) / s[mw,iw] * dt) * wave[:,mw,iw] ) - - if state.problem_data['arrival_state'] == True: - ## update q_hbox - for i in range(1, n_arrival_times): - q_hbox = q_hbox_initial.copy() - - for mw in range(num_waves): - if s[mw,iw-2] > 0: - q_hbox[:,0] -= arrival_times[i] / dx * (max(0, s[mw,iw-2] * arrival_times[i] - alpha * dx) / (s[mw,iw-2] * arrival_times[i]) * wave[:,mw,iw-2]) - - if s[mw, iw-1] < 0: - q_hbox[:,0] -= arrival_times[i] / dx * (min(-s[mw,iw-1] * arrival_times[i], (1 - alpha) * dx) / (-s[mw,iw-1] * arrival_times[i]) * wave[:,mw,iw-1]) - - if s_cells[mw,0,i] > 0: - for j in range(i): - q_hbox[:,0] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,0,j]) - - if s_cells[mw,0,i] * arrival_times[i] > alpha * dx: # - 1e-14: - # check the arrival wave - wave_cells[:,mw,0,i] = 0.0 - - if s_cells[mw,1,i] < 0: - for j in range(i): - q_hbox[:,0] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,1,j]) - - if s_cells[mw,1,i] > 0: - for j in range(i): - q_hbox[:,1] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,1,j]) - - if s_cells[mw,2,i] < 0: - for j in range(i): - q_hbox[:,1] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,2,j]) - - if (-s_cells[mw,2,i] * arrival_times[i]) > (1 - alpha) * dx: # - 1e-14: - # check the arrival wave - wave_cells[:,mw,2,i] = 0.0 - - if s[mw,iw+1] > 0: - q_hbox[:,1] -= arrival_times[i] / dx * (min(s[mw,iw+1] * arrival_times[i], alpha * dx) / (-s[mw,iw+1] * arrival_times[i]) * wave[:,mw,iw+1]) - - if s[mw,iw+2] < 0: - q_hbox[:,1] -= arrival_times[i] / dx * (max(0, -s[mw,iw+2] * arrival_times[i] - (1 - alpha) * dx) / (-s[mw,iw+2] * arrival_times[i]) * wave[:,mw,iw+2]) - - - - - wave_cells[:,:,1,i],s_cells[:,1,i],amdq_arr,apdq_arr = self.rp(q_hbox[:,0],q_hbox[:,1],aux_hbox[:,0],aux_hbox[:,1],state.problem_data) -#update wave cells :,:,0,i and others by doing middle value update every arrival step. make qhbox_middle -#add code in shallow_fwave_hbox_dry_1d with fw[iw-1] = middle and q_iw-1 and fw[iw+1] = middle and q_iw+2 - ## update q[iw-1], q[iw], q[iw+1] and q[iw+2] - arrival_times = np.append(arrival_times, dt) - n_arrival_times = len(arrival_times) - - for mw in range(num_waves): - for i in range(n_arrival_times-1): - if s_cells[mw,0,i] > 0: - q[:,iw] -= (arrival_times[i+1] - arrival_times[i]) / (alpha * dx) * (wave_cells[:,mw,0,i]) - if s_cells[mw,2,i] < 0: - q[:,iw+1] -= (arrival_times[i+1] - arrival_times[i]) / ((1 - alpha) * dx) * (wave_cells[:,mw,2,i]) - if s_cells[mw,1,i] < 0: - q[:,iw-1] -= (dt - arrival_times[i]) / dx * ( max(0, -s_cells[mw,1,i] * (dt - arrival_times[i]) - alpha * dx) / (-s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - q[:,iw] -= (dt - arrival_times[i]) / (alpha * dx) * ( min(-s_cells[mw,1, i] * (dt - arrival_times[i]), alpha * dx) / (-s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - if s_cells[mw,1,i] > 0: - q[:,iw+1] -= (dt - arrival_times[i]) / ((1 - alpha) * dx) * ( min(s_cells[mw,1, i] * (dt - arrival_times[i]), (1 - alpha) * dx) / (s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - q[:,iw+2] -= (dt - arrival_times[i]) / dx * ( max(0, s_cells[mw,1,i] * (dt - arrival_times[i]) - (1- alpha) * dx) / (s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - - - # Compute maximum wave speed - # add additional conditions for h-box - cfl = 0.0 - if 'method' not in state.problem_data: - for mw in range(wave.shape[1]): - smax1 = np.max(dtdx[LL:UL]*s[mw,LL-1:UL-1]) - smax2 = np.max(-dtdx[LL-1:UL-1]*s[mw,LL-1:UL-1]) - cfl = max(cfl,smax1,smax2) - elif state.problem_data['method'] == 'h_box': - for mw in range(wave.shape[1]): - smax1 = np.max(dtdx_hbox[LL:UL]*s[mw,LL-1:UL-1]) - smax2 = np.max(-dtdx_hbox[LL-1:UL-1]*s[mw,LL-1:UL-1]) - cfl = max(cfl,smax1,smax2) - elif state.problem_data['method'] == 'h_box_wave': - for mw in range(wave.shape[1]): - smax1 = np.max(dtdx_hbox[LL:UL]*s[mw,LL-1:UL-1]) - smax2 = np.max(-dtdx_hbox[LL-1:UL-1]*s[mw,LL-1:UL-1]) - cfl = max(cfl,smax1,smax2) - - - # If we are doing slope limiting we have more work to do - if self.order == 2: - # Initialize flux corrections - f = np.zeros( (num_eqn,grid.num_cells[0] + 2*self.num_ghost) ) - - # Apply Limiters to waves - if (limiter > 0).any(): - wave = tvd.limit(state.num_eqn,wave,s,limiter,dtdx) - - # Compute correction fluxes for second order q_{xx} terms - dtdxave = 0.5 * (dtdx[LL-1:UL-1] + dtdx[LL:UL]) - if self.fwave: - for mw in range(wave.shape[1]): - sabs = np.abs(s[mw,LL-1:UL-1]) - om = 1.0 - sabs*dtdxave[:UL-LL] - ssign = np.sign(s[mw,LL-1:UL-1]) - for m in range(num_eqn): - f[m,LL:UL] += 0.5 * ssign * om * wave[m,mw,LL-1:UL-1] - else: - for mw in range(wave.shape[1]): - sabs = np.abs(s[mw,LL-1:UL-1]) - om = 1.0 - sabs*dtdxave[:UL-LL] - for m in range(num_eqn): - f[m,LL:UL] += 0.5 * sabs * om * wave[m,mw,LL-1:UL-1] - - # Update q by differencing correction fluxes - for m in range(num_eqn): - q[m,LL:UL-1] -= dtdx[LL:UL-1] * (f[m,LL+1:UL] - f[m,LL:UL-1]) - - - else: raise Exception("Unrecognized kernel_language; choose 'Fortran' or 'Python'") - - self.cfl.update_global_max(cfl) - state.set_q_from_qbc(num_ghost,self.qbc) - if state.num_aux > 0: - state.set_aux_from_auxbc(num_ghost,self.auxbc) - - -# ============================================================================ -# ClawPack 2d Solver Class -# ============================================================================ -class ClawSolver2D(ClawSolver): - r""" - 2D Classic (Clawpack) solver. - - Solve using the wave propagation algorithms of Randy LeVeque's - Clawpack code (www.clawpack.org). - - In addition to the attributes of ClawSolver1D, ClawSolver2D - also has the following options: - - .. attribute:: dimensional_split - - If True, use dimensional splitting (Godunov splitting). - Dimensional splitting with Strang splitting is not supported - at present but could easily be enabled if necessary. - If False, use unsplit Clawpack algorithms, possibly including - transverse Riemann solves. - - .. attribute:: transverse_waves - - If dimensional_split is True, this option has no effect. If - dimensional_split is False, then transverse_waves should be one of - the following values: - - ClawSolver2D.no_trans: Transverse Riemann solver - not used. The stable CFL for this algorithm is 0.5. Not recommended. - - ClawSolver2D.trans_inc: Transverse increment waves are computed - and propagated. - - ClawSolver2D.trans_cor: Transverse increment waves and transverse - correction waves are computed and propagated. - - Note that only the fortran routines are supported for now in 2D. - """ - - __doc__ += add_parent_doc(ClawSolver) - - no_trans = 0 - trans_inc = 1 - trans_cor = 2 - - def __init__(self,riemann_solver=None, claw_package=None): - r""" - Create 2d Clawpack solver - - See :class:`ClawSolver2D` for more info. - """ - self.dimensional_split = True - self.transverse_waves = self.trans_inc - - self.num_dim = 2 - self.reflect_index = [1,2] - - self.aux1 = None - self.aux2 = None - self.aux3 = None - self.work = None - - super(ClawSolver2D,self).__init__(riemann_solver, claw_package) - - def _check_cfl_settings(self): - if (not self.dimensional_split) and (self.transverse_waves==0): - cfl_recommended = 0.5 - else: - cfl_recommended = 1.0 - - if self.cfl_max > cfl_recommended: - import warnings - warnings.warn('cfl_max is set higher than the recommended value of %s' % cfl_recommended) - warnings.warn(str(self.cfl_desired)) - - - def _allocate_workspace(self,solution): - r""" - Pack parameters into format recognized by Clawpack (Fortran) code. - - Sets the method array and the cparam common block for the Riemann solver. - """ - import numpy as np - - state = solution.state - - num_eqn,num_aux,num_waves,num_ghost,aux = state.num_eqn,state.num_aux,self.num_waves,self.num_ghost,state.aux - - #The following is a hack to work around an issue - #with f2py. It involves wastefully allocating three arrays. - #f2py seems not able to handle multiple zero-size arrays being passed. - # it appears the bug is related to f2py/src/fortranobject.c line 841. - if aux is None: num_aux=1 - - grid = state.grid - maxmx,maxmy = grid.num_cells[0],grid.num_cells[1] - maxm = max(maxmx, maxmy) - - # These work arrays really ought to live inside a fortran module - # as is done for sharpclaw - self.aux1 = np.empty((num_aux,maxm+2*num_ghost),order='F') - self.aux2 = np.empty((num_aux,maxm+2*num_ghost),order='F') - self.aux3 = np.empty((num_aux,maxm+2*num_ghost),order='F') - mwork = (maxm+2*num_ghost) * (5*num_eqn + num_waves + num_eqn*num_waves) - self.work = np.empty((mwork),order='F') - - - # ========== Hyperbolic Step ===================================== - def step_hyperbolic(self,solution): - r""" - Take a step on the homogeneous hyperbolic system using the Clawpack - algorithm. - - Clawpack is based on the Lax-Wendroff method, combined with Riemann - solvers and TVD limiters applied to waves. - """ - if(self.kernel_language == 'Fortran'): - state = solution.states[0] - grid = state.grid - dx,dy = grid.delta - mx,my = grid.num_cells - maxm = max(mx,my) - - self._apply_bcs(state) - qold = self.qbc.copy('F') - - rpn2 = self.rp.rpn2._cpointer - - if (self.dimensional_split) or (self.transverse_waves==0): - rpt2 = rpn2 # dummy value; it won't be called - else: - rpt2 = self.rp.rpt2._cpointer - - if self.dimensional_split: - #Right now only Godunov-dimensional-splitting is implemented. - #Strang-dimensional-splitting could be added following dimsp2.f in Clawpack. - - self.qbc, cfl_x = self.fmod.step2ds(maxm,self.num_ghost,mx,my, \ - qold,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,1,self.fwave,rpn2,rpt2) - - self.qbc, cfl_y = self.fmod.step2ds(maxm,self.num_ghost,mx,my, \ - self.qbc,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,2,self.fwave,rpn2,rpt2) - - cfl = max(cfl_x,cfl_y) - - else: - - self.qbc, cfl = self.fmod.step2(maxm,self.num_ghost,mx,my, \ - qold,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,self.fwave,rpn2,rpt2) - - self.cfl.update_global_max(cfl) - state.set_q_from_qbc(self.num_ghost,self.qbc) - if state.num_aux > 0: - state.set_aux_from_auxbc(self.num_ghost,self.auxbc) - - else: - raise NotImplementedError("No python implementation for step_hyperbolic in 2D.") - -# ============================================================================ -# ClawPack 3d Solver Class -# ============================================================================ -class ClawSolver3D(ClawSolver): - r""" - 3D Classic (Clawpack) solver. - - Solve using the wave propagation algorithms of Randy LeVeque's - Clawpack code (www.clawpack.org). - - In addition to the attributes of ClawSolver, ClawSolver3D - also has the following options: - - .. attribute:: dimensional_split - - If True, use dimensional splitting (Godunov splitting). - Dimensional splitting with Strang splitting is not supported - at present but could easily be enabled if necessary. - If False, use unsplit Clawpack algorithms, possibly including - transverse Riemann solves. - - .. attribute:: transverse_waves - - If dimensional_split is True, this option has no effect. If - dim_plit is False, then transverse_waves should be one of - the following values: - - ClawSolver3D.no_trans: Transverse Riemann solver - not used. The stable CFL for this algorithm is 0.5. Not recommended. - - ClawSolver3D.trans_inc: Transverse increment waves are computed - and propagated. - - ClawSolver3D.trans_cor: Transverse increment waves and transverse - correction waves are computed and propagated. - - Note that only Fortran routines are supported for now in 3D -- - there is no pure-python version. - """ - - __doc__ += add_parent_doc(ClawSolver) - - no_trans = 0 - trans_inc = 11 - trans_cor = 22 - - def __init__(self, riemann_solver=None, claw_package=None): - r""" - Create 3d Clawpack solver - - See :class:`ClawSolver3D` for more info. - """ - # Add the functions as required attributes - self.dimensional_split = True - self.transverse_waves = self.trans_cor - - self.num_dim = 3 - self.reflect_index = [1,2,3] - - self.aux1 = None - self.aux2 = None - self.aux3 = None - self.work = None - - super(ClawSolver3D,self).__init__(riemann_solver, claw_package) - - # ========== Setup routine ============================= - def _allocate_workspace(self,solution): - r""" - Allocate auxN and work arrays for use in Fortran subroutines. - """ - import numpy as np - - state = solution.states[0] - - num_eqn,num_aux,num_waves,num_ghost,aux = state.num_eqn,state.num_aux,self.num_waves,self.num_ghost,state.aux - - #The following is a hack to work around an issue - #with f2py. It involves wastefully allocating three arrays. - #f2py seems not able to handle multiple zero-size arrays being passed. - # it appears the bug is related to f2py/src/fortranobject.c line 841. - if(aux is None): num_aux=1 - - grid = state.grid - maxmx,maxmy,maxmz = grid.num_cells[0],grid.num_cells[1],grid.num_cells[2] - maxm = max(maxmx, maxmy, maxmz) - - # These work arrays really ought to live inside a fortran module - # as is done for sharpclaw - self.aux1 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') - self.aux2 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') - self.aux3 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') - mwork = (maxm+2*num_ghost) * (31*num_eqn + num_waves + num_eqn*num_waves) - self.work = np.empty((mwork),order='F') - - - # ========== Hyperbolic Step ===================================== - def step_hyperbolic(self,solution): - r""" - Take a step on the homogeneous hyperbolic system using the Clawpack - algorithm. - - Clawpack is based on the Lax-Wendroff method, combined with Riemann - solvers and TVD limiters applied to waves. - """ - if(self.kernel_language == 'Fortran'): - state = solution.states[0] - grid = state.grid - dx,dy,dz = grid.delta - mx,my,mz = grid.num_cells - maxm = max(mx,my,mz) - - self._apply_bcs(state) - qnew = self.qbc - qold = qnew.copy('F') - - rpn3 = self.rp.rpn3._cpointer - - if (self.dimensional_split) or (self.transverse_waves==0): - rpt3 = rpn3 # dummy value; it won't be called - rptt3 = rpn3 # dummy value; it won't be called - else: - rpt3 = self.rp.rpt3._cpointer - rptt3 = self.rp.rptt3._cpointer - - if self.dimensional_split: - #Right now only Godunov-dimensional-splitting is implemented. - #Strang-dimensional-splitting could be added following dimsp3.f in Clawpack. - - q, cfl_x = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ - qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,1,self.fwave,rpn3,rpt3,rptt3) - - q, cfl_y = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ - q,q,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,2,self.fwave,rpn3,rpt3,rptt3) - - q, cfl_z = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ - q,q,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,3,self.fwave,rpn3,rpt3,rptt3) - - cfl = max(cfl_x,cfl_y,cfl_z) - - else: - - q, cfl = self.fmod.step3(maxm,self.num_ghost,mx,my,mz, \ - qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,self.fwave,rpn3,rpt3,rptt3) - - self.cfl.update_global_max(cfl) - state.set_q_from_qbc(self.num_ghost,self.qbc) - if state.num_aux > 0: - state.set_aux_from_auxbc(self.num_ghost,self.auxbc) - - else: - raise NotImplementedError("No python implementation for step_hyperbolic in 3D.") From 7ef72b39b89bed3c678bc814a50003ec5f0c933e Mon Sep 17 00:00:00 2001 From: cr2940 Date: Fri, 18 Oct 2019 23:33:54 -0400 Subject: [PATCH 42/47] corrected test file for nonuniform advection 1d --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- examples/advection_1d/test_advection.py | 2 +- .../advection_1d/test_advection_nonunif.py | 24 +++++++++++-------- examples/shallow_1d/dam_break.py | 2 -- src/pyclaw/classic/solver.py | 4 ++-- 5 files changed, 19 insertions(+), 17 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 89567a0f8..120e0d0eb 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -33,8 +33,8 @@ def mapc2p_nonunif(xc): -def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='sharpclaw', weno_order=5, - time_integrator='SSPLMMk3', outdir='./_output'): +def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='sharpclaw', weno_order=17, + time_integrator='SSP104', outdir='./_output'): if use_petsc: import clawpack.petclaw as pyclaw diff --git a/examples/advection_1d/test_advection.py b/examples/advection_1d/test_advection.py index edd54bb8a..116041bfb 100644 --- a/examples/advection_1d/test_advection.py +++ b/examples/advection_1d/test_advection.py @@ -4,7 +4,7 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ - from . import advection_1d + import advection_1d def verify_expected(expected): """ given an expected value, returns a verification function """ diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index e33280b2a..0ee6db997 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -4,8 +4,8 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ - from . import advection_1d - from . import advection_1d_nonunif + import advection_1d + import advection_1d_nonunif def verify_expected(expected): """ given an expected value, returns a verification function """ @@ -42,7 +42,7 @@ def advection_nu_verify(claw): dx=claw.solution.domain.grid.delta[0] grid1d = claw.frames[0].state.grid grid1d.mapc2p = mapc2p_nonunif - nx = 500 + nx = 100 aux = np.zeros((1,nx)) aux[0,:] = np.diff(grid1d.p_nodes)/np.diff(grid1d.x.nodes) test = abs(np.sum(dx*aux[0,:]*(qfinal-q0))) @@ -53,25 +53,29 @@ def advection_nu_verify(claw): from clawpack.pyclaw.util import gen_variants - classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.156102048093389e-05), + classic_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.145595120502496e-16), kernel_languages=('Python','Fortran'), solver_type='classic', outdir=None) - sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.2345173603527554e-07), - kernel_languages=('Python','Fortran'), + sharp_tests_rk = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.677906041554036e-05), + kernel_languages=('Fortran',), # Not working for Python kernel solver_type='sharpclaw',time_integrator='SSP104', outdir=None) - sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.712477331998265e-07), - kernel_languages=('Python','Fortran'), + sharp_tests_lmm = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.6815055672658585e-08), + kernel_languages=('Python',), + solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) + + sharp_tests_lmm2 = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(1.6779069237771793e-05), + kernel_languages=('Fortran',), solver_type='sharpclaw',time_integrator='SSPLMMk3', outdir=None) - weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(4.014894137206369e-08), + weno_tests = gen_variants(advection_1d_nonunif.setup, verify_expected_nonunif(2.4413818316872576e-04), kernel_languages=('Fortran',), solver_type='sharpclaw', time_integrator='SSP104', weno_order=17, outdir=None) from itertools import chain - for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests): + for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, sharp_tests_lmm2, weno_tests): yield test diff --git a/examples/shallow_1d/dam_break.py b/examples/shallow_1d/dam_break.py index cb09ad54a..ff79cd18a 100755 --- a/examples/shallow_1d/dam_break.py +++ b/examples/shallow_1d/dam_break.py @@ -87,8 +87,6 @@ def setup(use_petsc=False,kernel_language='Fortran',outdir='./_output',solver_ty eps=0.1 state.q[depth,:] = 1.0 + eps*np.exp(-(xc-x0)**2/0.5) state.q[momentum,:] = 0. - state.aux[0,:] = np.zeros((1,mx)) - state.aux[0,int(mx/2):] +=0.5 claw = pyclaw.Controller() claw.keep_copy = True diff --git a/src/pyclaw/classic/solver.py b/src/pyclaw/classic/solver.py index c4a6c2b4d..401a02620 100644 --- a/src/pyclaw/classic/solver.py +++ b/src/pyclaw/classic/solver.py @@ -79,7 +79,7 @@ def __init__(self,riemann_solver=None,claw_package=None): """ self.num_ghost = 2 self.limiters = tvd.minmod - self.order = 1 + self.order = 2 self.source_split = 1 self.fwave = False self.step_source = None @@ -289,7 +289,6 @@ def step_hyperbolic(self,solution): state = solution.states[0] grid = state.grid - #print(grid) self._apply_bcs(state) @@ -307,6 +306,7 @@ def step_hyperbolic(self,solution): q = self.qbc aux = self.auxbc + # Limiter to use in the pth family limiter = np.array(self._mthlim,ndmin=1) dtdx = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) From a5aafd30b6570421772d41d2fe16fc014393c1d9 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Fri, 18 Oct 2019 23:40:05 -0400 Subject: [PATCH 43/47] __init__.py for advection 1d --- examples/advection_1d/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 examples/advection_1d/__init__.py diff --git a/examples/advection_1d/__init__.py b/examples/advection_1d/__init__.py new file mode 100644 index 000000000..e69de29bb From 7025949445191ed7a70d525557d1de666a19e3ae Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Fri, 18 Oct 2019 23:48:20 -0400 Subject: [PATCH 44/47] Update advection_1d_nonunif.py set solver to classic and kernel to python --- examples/advection_1d/advection_1d_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 120e0d0eb..b156b46fd 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -33,7 +33,7 @@ def mapc2p_nonunif(xc): -def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='sharpclaw', weno_order=17, +def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output'): if use_petsc: From 2ca11cc6e64e71a27c761a7ab201061b2be55bbd Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Fri, 18 Oct 2019 23:48:47 -0400 Subject: [PATCH 45/47] Update advection_1d_nonunif.py --- examples/advection_1d/advection_1d_nonunif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index b156b46fd..b9b33f921 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -33,7 +33,7 @@ def mapc2p_nonunif(xc): -def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='classic', weno_order=5, +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output'): if use_petsc: From 97aa84b8f882298ecbf7bde5c9f95dbb9f49aa91 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Mon, 21 Oct 2019 17:21:42 -0400 Subject: [PATCH 46/47] get rid of plots and output directories --- examples/advection_1d/advection_1d_nonunif.py | 4 ++-- examples/advection_1d/test_advection.py | 2 +- examples/advection_1d/test_advection_nonunif.py | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 120e0d0eb..c40c89ea3 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -33,7 +33,7 @@ def mapc2p_nonunif(xc): -def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='sharpclaw', weno_order=17, +def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output'): if use_petsc: @@ -112,7 +112,7 @@ def setplot(plotdata): # Set up for axes in this figure: plotaxes = plotfigure.new_plotaxes() plotaxes.xlimits = [-0.25,0.25] - plotaxes.ylimits = [-.2,1.0] + plotaxes.ylimits = [-.2,1.5] plotaxes.title = 'q' # Set up for item on these axes: diff --git a/examples/advection_1d/test_advection.py b/examples/advection_1d/test_advection.py index 116041bfb..edd54bb8a 100644 --- a/examples/advection_1d/test_advection.py +++ b/examples/advection_1d/test_advection.py @@ -4,7 +4,7 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ - import advection_1d + from . import advection_1d def verify_expected(expected): """ given an expected value, returns a verification function """ diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py index 0ee6db997..d7c0df411 100644 --- a/examples/advection_1d/test_advection_nonunif.py +++ b/examples/advection_1d/test_advection_nonunif.py @@ -4,8 +4,8 @@ def test_1d_advection(): tests against expected classic, sharpclaw, and high-order weno results """ - import advection_1d - import advection_1d_nonunif + from . import advection_1d + from . import advection_1d_nonunif def verify_expected(expected): """ given an expected value, returns a verification function """ From 98d2549d948c773c6188e820885c94c52b6e41ef Mon Sep 17 00:00:00 2001 From: cr2940 Date: Tue, 22 Oct 2019 12:38:29 -0400 Subject: [PATCH 47/47] clean up examples and added nonuniform advection 1d example --- examples/advection_1d/__init__.py | 0 src/pyclaw/classic/solverjli.py | 866 ------------------------------ 2 files changed, 866 deletions(-) create mode 100644 examples/advection_1d/__init__.py delete mode 100644 src/pyclaw/classic/solverjli.py diff --git a/examples/advection_1d/__init__.py b/examples/advection_1d/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/pyclaw/classic/solverjli.py b/src/pyclaw/classic/solverjli.py deleted file mode 100644 index 7dfd80c57..000000000 --- a/src/pyclaw/classic/solverjli.py +++ /dev/null @@ -1,866 +0,0 @@ -r""" -Module containing the classic Clawpack solvers. - -This module contains the pure and wrapped classic clawpack solvers. All -clawpack solvers inherit from the :class:`ClawSolver` superclass which in turn -inherits from the :class:`~pyclaw.solver.Solver` superclass. These -are both pure virtual classes; the only solver classes that should be instantiated -are the dimension-specific ones, :class:`ClawSolver1D` and :class:`ClawSolver2D`. -""" - -from clawpack.pyclaw.util import add_parent_doc -from clawpack.pyclaw.solver import Solver -from clawpack.pyclaw.limiters import tvd - -# ============================================================================ -# Generic Clawpack solver class -# ============================================================================ -class ClawSolver(Solver): - r""" - Generic classic Clawpack solver - - All Clawpack solvers inherit from this base class. - - .. attribute:: mthlim - - Limiter(s) to be used. Specified either as one value or a list. - If one value, the specified limiter is used for all wave families. - If a list, the specified values indicate which limiter to apply to - each wave family. Take a look at pyclaw.limiters.tvd for an enumeration. - ``Default = limiters.tvd.minmod`` - - .. attribute:: order - - Order of the solver, either 1 for first order (i.e., Godunov's method) - or 2 for second order (Lax-Wendroff-LeVeque). - ``Default = 2`` - - .. attribute:: source_split - - Which source splitting method to use: 1 for first - order Godunov splitting and 2 for second order Strang splitting. - ``Default = 1`` - - .. attribute:: fwave - - Whether to split the flux jump (rather than the jump in Q) into waves; - requires that the Riemann solver performs the splitting. - ``Default = False`` - - .. attribute:: step_source - - Handle for function that evaluates the source term. - The required signature for this function is: - - def step_source(solver,state,dt) - - .. attribute:: kernel_language - - Specifies whether to use wrapped Fortran routines ('Fortran') - or pure Python ('Python'). ``Default = 'Fortran'``. - - .. attribute:: verbosity - - The level of detail of logged messages from the Fortran solver. - ``Default = 0``. - - """ - - # ========== Generic Init Routine ======================================== - def __init__(self,riemann_solver=None,claw_package=None): - r""" - See :class:`ClawSolver` for full documentation. - - Output: - - (:class:`ClawSolver`) - Initialized clawpack solver - """ - self.num_ghost = 2 - self.limiters = tvd.vanleer - self.order = 2 - self.source_split = 1 - self.fwave = False - self.step_source = None - self.kernel_language = 'Fortran' - self.verbosity = 0 - self.cfl_max = 1.0 - self.cfl_desired = 0.9 - self._mthlim = self.limiters - self._method = None - self.dt_old = None - - # Call general initialization function - super(ClawSolver,self).__init__(riemann_solver,claw_package) - - # ========== Time stepping routines ====================================== - def step(self,solution,take_one_step,tstart,tend): - r""" - Evolve solution one time step - - The elements of the algorithm for taking one step are: - - 1. Pick a step size as specified by the base solver attribute :func:`get_dt` - - 2. A half step on the source term :func:`step_source` if Strang splitting is - being used (:attr:`source_split` = 2) - - 3. A step on the homogeneous problem :math:`q_t + f(q)_x = 0` is taken - - 4. A second half step or a full step is taken on the source term - :func:`step_source` depending on whether Strang splitting was used - (:attr:`source_split` = 2) or Godunov splitting - (:attr:`source_split` = 1) - - This routine is called from the method evolve_to_time defined in the - pyclaw.solver.Solver superclass. - - :Input: - - *solution* - (:class:`~pyclaw.solution.Solution`) solution to be evolved - - :Output: - - (bool) - True if full step succeeded, False otherwise - """ - self.get_dt(solution.t,tstart,tend,take_one_step) - self.cfl.set_global_max(0.) - - if self.source_split == 2 and self.step_source is not None: - self.step_source(self,solution.states[0],self.dt/2.0) - - self.step_hyperbolic(solution) - - # Check here if the CFL condition is satisfied. - # If not, return # immediately to evolve_to_time and let it deal with - # picking a new step size (dt). - if self.cfl.get_cached_max() >= self.cfl_max: - return False - - if self.step_source is not None: - # Strang splitting - if self.source_split == 2: - self.step_source(self,solution.states[0],self.dt/2.0) - - # Godunov Splitting - if self.source_split == 1: - self.step_source(self,solution.states[0],self.dt) - - return True - - def _check_cfl_settings(self): - pass - - def _allocate_workspace(self,solution): - pass - - def step_hyperbolic(self,solution): - r""" - Take one homogeneous step on the solution. - - This is a dummy routine and must be overridden. - """ - raise Exception("Dummy routine, please override!") - - def _set_mthlim(self): - r""" - Convenience routine to convert users limiter specification to - the format understood by the Fortran code (i.e., a list of length num_waves). - """ - #self._mthlim = self.limiters - #if not isinstance(self.limiters,list): self._mthlim=[self._mthlim] - # if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves - # if len(self._mthlim)!=self.num_waves: - # raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves') - - def _set_method(self,state): - r""" - Set values of the solver._method array required by the Fortran code. - These are algorithmic parameters. - """ - import numpy as np - #We ought to put method and many other things in a Fortran - #module and set the fortran variables directly here. - self._method =np.empty(7, dtype=int,order='F') - self._method[0] = self.dt_variable - self._method[1] = self.order - if self.num_dim==1: - self._method[2] = 0 # Not used in 1D - elif self.dimensional_split: - self._method[2] = -1 # First-order dimensional splitting - else: - self._method[2] = self.transverse_waves - self._method[3] = self.verbosity - self._method[4] = 0 # Not used for PyClaw (would be self.source_split) - self._method[5] = state.index_capa + 1 - self._method[6] = state.num_aux - - def setup(self,solution): - r""" - Perform essential solver setup. This routine must be called before - solver.step() may be called. - """ - # This is a hack to deal with the fact that petsc4py - # doesn't allow us to change the stencil_width (num_ghost) - solution.state.set_num_ghost(self.num_ghost) - # End hack - - self._check_cfl_settings() - - self._set_mthlim() - if(self.kernel_language == 'Fortran'): - if self.fmod is None: - so_name = 'clawpack.pyclaw.classic.classic'+str(self.num_dim) - self.fmod = __import__(so_name,fromlist=['clawpack.pyclaw.classic']) - self._set_fortran_parameters(solution) - self._allocate_workspace(solution) - elif self.num_dim>1: - raise Exception('Only Fortran kernels are supported in multi-D.') - - self._allocate_bc_arrays(solution.states[0]) - - super(ClawSolver,self).setup(solution) - - - def _set_fortran_parameters(self,solution): - r""" - Pack parameters into format recognized by Clawpack (Fortran) code. - - Sets the solver._method array and the cparam common block for the Riemann solver. - """ - self._set_method(solution.state) - # The reload here is necessary because otherwise the common block - # cparam in the Riemann solver doesn't get flushed between running - # different tests in a single Python session. - reload(self.fmod) - solution.state.set_cparam(self.fmod) - solution.state.set_cparam(self.rp) - - def __del__(self): - r""" - Delete Fortran objects, which otherwise tend to persist in Python sessions. - """ - if(self.kernel_language == 'Fortran'): - del self.fmod - - super(ClawSolver,self).__del__() - - -# ============================================================================ -# ClawPack 1d Solver Class -# ============================================================================ -class ClawSolver1D(ClawSolver): - r""" - Clawpack evolution routine in 1D - - This class represents the 1d clawpack solver on a single grid. Note that - there are routines here for interfacing with the fortran time stepping - routines and the Python time stepping routines. The ones used are - dependent on the argument given to the initialization of the solver - (defaults to python). - - """ - - __doc__ += add_parent_doc(ClawSolver) - - def __init__(self, riemann_solver=None, claw_package=None): - r""" - Create 1d Clawpack solver - - Output: - - (:class:`ClawSolver1D`) - Initialized 1d clawpack solver - - See :class:`ClawSolver1D` for more info. - """ - self.num_dim = 1 - self.reflect_index = [1] - - super(ClawSolver1D,self).__init__(riemann_solver, claw_package) - - - # ========== Homogeneous Step ===================================== - def step_hyperbolic(self,solution): - r""" - Take one time step on the homogeneous hyperbolic system. - - :Input: - - *solution* - (:class:`~pyclaw.solution.Solution`) Solution that - will be evolved - """ - import numpy as np - - state = solution.states[0] - grid = state.grid - - - self._apply_bcs(state) - - num_eqn,num_ghost = state.num_eqn,self.num_ghost - - if(self.kernel_language == 'Fortran'): - mx = grid.num_cells[0] - dx,dt = grid.delta[0],self.dt - dtdx = np.zeros( (mx+2*num_ghost) ) + dt/dx - rp1 = self.rp.rp1._cpointer - - self.qbc,cfl = self.fmod.step1(num_ghost,mx,self.qbc,self.auxbc,dx,dt,self._method,self._mthlim,self.fwave,rp1) - - elif(self.kernel_language == 'Python'): - q = self.qbc - aux = self.auxbc - # Limiter to use in the pth family - limiter = np.array(self._mthlim,ndmin=1) - - dtdx = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - - # Find local value for dt/dx - if 'method' not in state.problem_data: - if state.index_capa>=0: - dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) - else: - dtdx += self.dt / grid.delta[0] - elif state.problem_data['method'] == 'h_box': - xpxc = state.problem_data['xpxc'] - dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) - dtdx_hbox = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - dtdx_hbox += self.dt / (grid.delta[0] * xpxc) - elif state.problem_data['method'] == 'h_box_wave': - xpxc = state.problem_data['xpxc'] - dtdx = self.dt / (grid.delta[0] * aux[state.index_capa,:]) - dtdx_hbox = np.zeros( (2*self.num_ghost+grid.num_cells[0]) ) - dtdx_hbox += self.dt / (grid.delta[0] * xpxc) - - - - # Solve Riemann problem at each interface - q_l=q[:,:-1].copy() - q_r=q[:,1:].copy() - if state.aux is not None: - aux_l=aux[:,:-1].copy() - aux_r=aux[:,1:].copy() - else: - aux_l = None - aux_r = None - - if 'method' not in state.problem_data: - # normal case - wave,s,amdq,apdq = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - - - # Update loop limits, these are the limits for the Riemann solver - # locations, which then update a grid cell value - # We include the Riemann problem just outside of the grid so we can - # do proper limiting at the grid edges - # LL | | UL - # | LL | | | | ... | | | UL | | - # | | - - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - - # Update q for Godunov update - for m in range(num_eqn): - q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] - q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] - - elif state.problem_data['method'] == 'h_box': - # # add corrections - wave,s,amdq,apdq,f_corr_l,f_corr_r = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - - # Update q for Godunov update - for m in range(num_eqn): - q[m,LL:UL] -= dtdx[LL:UL]*(apdq[m,LL-1:UL-1] - f_corr_r[m,LL-1:UL-1]) - q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*(amdq[m,LL-1:UL-1] + f_corr_l[m,LL-1:UL-1]) - - elif state.problem_data['method'] == 'h_box_wave': - # # add corrections - state.problem_data['arrival_state'] = False - - - wave,s,amdq,apdq,q_hbox_initial,aux_hbox = self.rp(q_l,q_r,aux_l,aux_r,state.problem_data) - LL = self.num_ghost - 1 - UL = self.num_ghost + grid.num_cells[0] + 1 - - # Update q for Godunov update - iw = state.problem_data['wall_position'] + self.num_ghost - 1 - q_last = q[:,iw:iw+2].copy() - - for m in range(num_eqn): - q[m,LL:UL] -= dtdx[LL:UL]*apdq[m,LL-1:UL-1] - q[m,LL-1:UL-1] -= dtdx[LL-1:UL-1]*amdq[m,LL-1:UL-1] - - # check the arrivals - q[:,iw:iw+2] = q_last[:,:] # reset the wall cells - dt = self.dt - num_waves = self.num_waves - dx = grid.delta[0] * xpxc - alpha = state.problem_data['fraction'] - arrival_times = np.array([0.0]) - for mw in range(num_waves): - if (s[mw,iw-1] > 0 and (s[mw,iw-1] * dt > alpha * dx)): - arrival_times = np.append(arrival_times, alpha * dx / s[mw,iw-1]) - if (s[mw,iw+1] < 0 and ( (-s[mw,iw+1]) * dt > (1 - alpha) * dx ) ): - arrival_times = np.append(arrival_times, -(1 - alpha) * dx / s[mw,iw+1]) - arrival_times.sort() - - n_arrival_times = len(arrival_times) - if n_arrival_times == 1 : - state.problem_data['arrival_state'] = False - else: - state.problem_data['arrival_state'] = True - s_cells = np.zeros((num_waves, 3, n_arrival_times)) - s_cells[:,:,0] = s[:, iw-1:iw+2].copy() - wave_cells = np.zeros((num_eqn, num_waves, 3, n_arrival_times)) - wave_cells[:,:,:,0] = wave[:,:,iw-1:iw+2].copy() - - - if state.problem_data['arrival_state'] == False: - q[:,iw] -= dt/(alpha * dx) * apdq[:,iw-1] - q[:,iw+1] -= dt/((1 - alpha)*dx) * amdq[:,iw+1] - for mw in range(num_waves): - if (s[mw,iw] < 0): - q[:,iw-1] -= dt/dx * ( max(0, -s[mw,iw] * dt - alpha * dx) / (-s[mw,iw] * dt) * wave[:,mw,iw] ) - q[:,iw] -= dt/(alpha * dx) * (min(-s[mw,iw] * dt, alpha * dx) / (-s[mw,iw] * dt) * wave[:,mw,iw] ) - elif (s[mw,iw] > 0): - q[:,iw+1] -= dt/((1 - alpha)*dx) * (min(s[mw,iw] * dt, (1 - alpha) * dx) / (s[mw,iw] * dt) * wave[:,mw,iw] ) - q[:,iw+2] -= dt/dx * ( (max(0, s[mw,iw] * dt - (1 - alpha) * dx) / s[mw,iw] * dt) * wave[:,mw,iw] ) - - if state.problem_data['arrival_state'] == True: - ## update q_hbox - for i in range(1, n_arrival_times): - q_hbox = q_hbox_initial.copy() - - for mw in range(num_waves): - if s[mw,iw-2] > 0: - q_hbox[:,0] -= arrival_times[i] / dx * (max(0, s[mw,iw-2] * arrival_times[i] - alpha * dx) / (s[mw,iw-2] * arrival_times[i]) * wave[:,mw,iw-2]) - - if s[mw, iw-1] < 0: - q_hbox[:,0] -= arrival_times[i] / dx * (min(-s[mw,iw-1] * arrival_times[i], (1 - alpha) * dx) / (-s[mw,iw-1] * arrival_times[i]) * wave[:,mw,iw-1]) - - if s_cells[mw,0,i] > 0: - for j in range(i): - q_hbox[:,0] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,0,j]) - - if s_cells[mw,0,i] * arrival_times[i] > alpha * dx: # - 1e-14: - # check the arrival wave - wave_cells[:,mw,0,i] = 0.0 - - if s_cells[mw,1,i] < 0: - for j in range(i): - q_hbox[:,0] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,1,j]) - - if s_cells[mw,1,i] > 0: - for j in range(i): - q_hbox[:,1] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,1,j]) - - if s_cells[mw,2,i] < 0: - for j in range(i): - q_hbox[:,1] -= (arrival_times[j+1] - arrival_times[j]) / dx * (wave_cells[:,mw,2,j]) - - if (-s_cells[mw,2,i] * arrival_times[i]) > (1 - alpha) * dx: # - 1e-14: - # check the arrival wave - wave_cells[:,mw,2,i] = 0.0 - - if s[mw,iw+1] > 0: - q_hbox[:,1] -= arrival_times[i] / dx * (min(s[mw,iw+1] * arrival_times[i], alpha * dx) / (-s[mw,iw+1] * arrival_times[i]) * wave[:,mw,iw+1]) - - if s[mw,iw+2] < 0: - q_hbox[:,1] -= arrival_times[i] / dx * (max(0, -s[mw,iw+2] * arrival_times[i] - (1 - alpha) * dx) / (-s[mw,iw+2] * arrival_times[i]) * wave[:,mw,iw+2]) - - - - - wave_cells[:,:,1,i],s_cells[:,1,i],amdq_arr,apdq_arr = self.rp(q_hbox[:,0],q_hbox[:,1],aux_hbox[:,0],aux_hbox[:,1],state.problem_data) -#update wave cells :,:,0,i and others by doing middle value update every arrival step. make qhbox_middle -#add code in shallow_fwave_hbox_dry_1d with fw[iw-1] = middle and q_iw-1 and fw[iw+1] = middle and q_iw+2 - ## update q[iw-1], q[iw], q[iw+1] and q[iw+2] - arrival_times = np.append(arrival_times, dt) - n_arrival_times = len(arrival_times) - - for mw in range(num_waves): - for i in range(n_arrival_times-1): - if s_cells[mw,0,i] > 0: - q[:,iw] -= (arrival_times[i+1] - arrival_times[i]) / (alpha * dx) * (wave_cells[:,mw,0,i]) - if s_cells[mw,2,i] < 0: - q[:,iw+1] -= (arrival_times[i+1] - arrival_times[i]) / ((1 - alpha) * dx) * (wave_cells[:,mw,2,i]) - if s_cells[mw,1,i] < 0: - q[:,iw-1] -= (dt - arrival_times[i]) / dx * ( max(0, -s_cells[mw,1,i] * (dt - arrival_times[i]) - alpha * dx) / (-s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - q[:,iw] -= (dt - arrival_times[i]) / (alpha * dx) * ( min(-s_cells[mw,1, i] * (dt - arrival_times[i]), alpha * dx) / (-s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - if s_cells[mw,1,i] > 0: - q[:,iw+1] -= (dt - arrival_times[i]) / ((1 - alpha) * dx) * ( min(s_cells[mw,1, i] * (dt - arrival_times[i]), (1 - alpha) * dx) / (s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - q[:,iw+2] -= (dt - arrival_times[i]) / dx * ( max(0, s_cells[mw,1,i] * (dt - arrival_times[i]) - (1- alpha) * dx) / (s_cells[mw,1,i] * (dt - arrival_times[i])) * wave_cells[:,mw,1,i] ) - - - # Compute maximum wave speed - # add additional conditions for h-box - cfl = 0.0 - if 'method' not in state.problem_data: - for mw in range(wave.shape[1]): - smax1 = np.max(dtdx[LL:UL]*s[mw,LL-1:UL-1]) - smax2 = np.max(-dtdx[LL-1:UL-1]*s[mw,LL-1:UL-1]) - cfl = max(cfl,smax1,smax2) - elif state.problem_data['method'] == 'h_box': - for mw in range(wave.shape[1]): - smax1 = np.max(dtdx_hbox[LL:UL]*s[mw,LL-1:UL-1]) - smax2 = np.max(-dtdx_hbox[LL-1:UL-1]*s[mw,LL-1:UL-1]) - cfl = max(cfl,smax1,smax2) - elif state.problem_data['method'] == 'h_box_wave': - for mw in range(wave.shape[1]): - smax1 = np.max(dtdx_hbox[LL:UL]*s[mw,LL-1:UL-1]) - smax2 = np.max(-dtdx_hbox[LL-1:UL-1]*s[mw,LL-1:UL-1]) - cfl = max(cfl,smax1,smax2) - - - # If we are doing slope limiting we have more work to do - if self.order == 2: - # Initialize flux corrections - f = np.zeros( (num_eqn,grid.num_cells[0] + 2*self.num_ghost) ) - - # Apply Limiters to waves - if (limiter > 0).any(): - wave = tvd.limit(state.num_eqn,wave,s,limiter,dtdx) - - # Compute correction fluxes for second order q_{xx} terms - dtdxave = 0.5 * (dtdx[LL-1:UL-1] + dtdx[LL:UL]) - if self.fwave: - for mw in range(wave.shape[1]): - sabs = np.abs(s[mw,LL-1:UL-1]) - om = 1.0 - sabs*dtdxave[:UL-LL] - ssign = np.sign(s[mw,LL-1:UL-1]) - for m in range(num_eqn): - f[m,LL:UL] += 0.5 * ssign * om * wave[m,mw,LL-1:UL-1] - else: - for mw in range(wave.shape[1]): - sabs = np.abs(s[mw,LL-1:UL-1]) - om = 1.0 - sabs*dtdxave[:UL-LL] - for m in range(num_eqn): - f[m,LL:UL] += 0.5 * sabs * om * wave[m,mw,LL-1:UL-1] - - # Update q by differencing correction fluxes - for m in range(num_eqn): - q[m,LL:UL-1] -= dtdx[LL:UL-1] * (f[m,LL+1:UL] - f[m,LL:UL-1]) - - - else: raise Exception("Unrecognized kernel_language; choose 'Fortran' or 'Python'") - - self.cfl.update_global_max(cfl) - state.set_q_from_qbc(num_ghost,self.qbc) - if state.num_aux > 0: - state.set_aux_from_auxbc(num_ghost,self.auxbc) - - -# ============================================================================ -# ClawPack 2d Solver Class -# ============================================================================ -class ClawSolver2D(ClawSolver): - r""" - 2D Classic (Clawpack) solver. - - Solve using the wave propagation algorithms of Randy LeVeque's - Clawpack code (www.clawpack.org). - - In addition to the attributes of ClawSolver1D, ClawSolver2D - also has the following options: - - .. attribute:: dimensional_split - - If True, use dimensional splitting (Godunov splitting). - Dimensional splitting with Strang splitting is not supported - at present but could easily be enabled if necessary. - If False, use unsplit Clawpack algorithms, possibly including - transverse Riemann solves. - - .. attribute:: transverse_waves - - If dimensional_split is True, this option has no effect. If - dimensional_split is False, then transverse_waves should be one of - the following values: - - ClawSolver2D.no_trans: Transverse Riemann solver - not used. The stable CFL for this algorithm is 0.5. Not recommended. - - ClawSolver2D.trans_inc: Transverse increment waves are computed - and propagated. - - ClawSolver2D.trans_cor: Transverse increment waves and transverse - correction waves are computed and propagated. - - Note that only the fortran routines are supported for now in 2D. - """ - - __doc__ += add_parent_doc(ClawSolver) - - no_trans = 0 - trans_inc = 1 - trans_cor = 2 - - def __init__(self,riemann_solver=None, claw_package=None): - r""" - Create 2d Clawpack solver - - See :class:`ClawSolver2D` for more info. - """ - self.dimensional_split = True - self.transverse_waves = self.trans_inc - - self.num_dim = 2 - self.reflect_index = [1,2] - - self.aux1 = None - self.aux2 = None - self.aux3 = None - self.work = None - - super(ClawSolver2D,self).__init__(riemann_solver, claw_package) - - def _check_cfl_settings(self): - if (not self.dimensional_split) and (self.transverse_waves==0): - cfl_recommended = 0.5 - else: - cfl_recommended = 1.0 - - if self.cfl_max > cfl_recommended: - import warnings - warnings.warn('cfl_max is set higher than the recommended value of %s' % cfl_recommended) - warnings.warn(str(self.cfl_desired)) - - - def _allocate_workspace(self,solution): - r""" - Pack parameters into format recognized by Clawpack (Fortran) code. - - Sets the method array and the cparam common block for the Riemann solver. - """ - import numpy as np - - state = solution.state - - num_eqn,num_aux,num_waves,num_ghost,aux = state.num_eqn,state.num_aux,self.num_waves,self.num_ghost,state.aux - - #The following is a hack to work around an issue - #with f2py. It involves wastefully allocating three arrays. - #f2py seems not able to handle multiple zero-size arrays being passed. - # it appears the bug is related to f2py/src/fortranobject.c line 841. - if aux is None: num_aux=1 - - grid = state.grid - maxmx,maxmy = grid.num_cells[0],grid.num_cells[1] - maxm = max(maxmx, maxmy) - - # These work arrays really ought to live inside a fortran module - # as is done for sharpclaw - self.aux1 = np.empty((num_aux,maxm+2*num_ghost),order='F') - self.aux2 = np.empty((num_aux,maxm+2*num_ghost),order='F') - self.aux3 = np.empty((num_aux,maxm+2*num_ghost),order='F') - mwork = (maxm+2*num_ghost) * (5*num_eqn + num_waves + num_eqn*num_waves) - self.work = np.empty((mwork),order='F') - - - # ========== Hyperbolic Step ===================================== - def step_hyperbolic(self,solution): - r""" - Take a step on the homogeneous hyperbolic system using the Clawpack - algorithm. - - Clawpack is based on the Lax-Wendroff method, combined with Riemann - solvers and TVD limiters applied to waves. - """ - if(self.kernel_language == 'Fortran'): - state = solution.states[0] - grid = state.grid - dx,dy = grid.delta - mx,my = grid.num_cells - maxm = max(mx,my) - - self._apply_bcs(state) - qold = self.qbc.copy('F') - - rpn2 = self.rp.rpn2._cpointer - - if (self.dimensional_split) or (self.transverse_waves==0): - rpt2 = rpn2 # dummy value; it won't be called - else: - rpt2 = self.rp.rpt2._cpointer - - if self.dimensional_split: - #Right now only Godunov-dimensional-splitting is implemented. - #Strang-dimensional-splitting could be added following dimsp2.f in Clawpack. - - self.qbc, cfl_x = self.fmod.step2ds(maxm,self.num_ghost,mx,my, \ - qold,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,1,self.fwave,rpn2,rpt2) - - self.qbc, cfl_y = self.fmod.step2ds(maxm,self.num_ghost,mx,my, \ - self.qbc,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,2,self.fwave,rpn2,rpt2) - - cfl = max(cfl_x,cfl_y) - - else: - - self.qbc, cfl = self.fmod.step2(maxm,self.num_ghost,mx,my, \ - qold,self.qbc,self.auxbc,dx,dy,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,self.fwave,rpn2,rpt2) - - self.cfl.update_global_max(cfl) - state.set_q_from_qbc(self.num_ghost,self.qbc) - if state.num_aux > 0: - state.set_aux_from_auxbc(self.num_ghost,self.auxbc) - - else: - raise NotImplementedError("No python implementation for step_hyperbolic in 2D.") - -# ============================================================================ -# ClawPack 3d Solver Class -# ============================================================================ -class ClawSolver3D(ClawSolver): - r""" - 3D Classic (Clawpack) solver. - - Solve using the wave propagation algorithms of Randy LeVeque's - Clawpack code (www.clawpack.org). - - In addition to the attributes of ClawSolver, ClawSolver3D - also has the following options: - - .. attribute:: dimensional_split - - If True, use dimensional splitting (Godunov splitting). - Dimensional splitting with Strang splitting is not supported - at present but could easily be enabled if necessary. - If False, use unsplit Clawpack algorithms, possibly including - transverse Riemann solves. - - .. attribute:: transverse_waves - - If dimensional_split is True, this option has no effect. If - dim_plit is False, then transverse_waves should be one of - the following values: - - ClawSolver3D.no_trans: Transverse Riemann solver - not used. The stable CFL for this algorithm is 0.5. Not recommended. - - ClawSolver3D.trans_inc: Transverse increment waves are computed - and propagated. - - ClawSolver3D.trans_cor: Transverse increment waves and transverse - correction waves are computed and propagated. - - Note that only Fortran routines are supported for now in 3D -- - there is no pure-python version. - """ - - __doc__ += add_parent_doc(ClawSolver) - - no_trans = 0 - trans_inc = 11 - trans_cor = 22 - - def __init__(self, riemann_solver=None, claw_package=None): - r""" - Create 3d Clawpack solver - - See :class:`ClawSolver3D` for more info. - """ - # Add the functions as required attributes - self.dimensional_split = True - self.transverse_waves = self.trans_cor - - self.num_dim = 3 - self.reflect_index = [1,2,3] - - self.aux1 = None - self.aux2 = None - self.aux3 = None - self.work = None - - super(ClawSolver3D,self).__init__(riemann_solver, claw_package) - - # ========== Setup routine ============================= - def _allocate_workspace(self,solution): - r""" - Allocate auxN and work arrays for use in Fortran subroutines. - """ - import numpy as np - - state = solution.states[0] - - num_eqn,num_aux,num_waves,num_ghost,aux = state.num_eqn,state.num_aux,self.num_waves,self.num_ghost,state.aux - - #The following is a hack to work around an issue - #with f2py. It involves wastefully allocating three arrays. - #f2py seems not able to handle multiple zero-size arrays being passed. - # it appears the bug is related to f2py/src/fortranobject.c line 841. - if(aux is None): num_aux=1 - - grid = state.grid - maxmx,maxmy,maxmz = grid.num_cells[0],grid.num_cells[1],grid.num_cells[2] - maxm = max(maxmx, maxmy, maxmz) - - # These work arrays really ought to live inside a fortran module - # as is done for sharpclaw - self.aux1 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') - self.aux2 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') - self.aux3 = np.empty((num_aux,maxm+2*num_ghost,3),order='F') - mwork = (maxm+2*num_ghost) * (31*num_eqn + num_waves + num_eqn*num_waves) - self.work = np.empty((mwork),order='F') - - - # ========== Hyperbolic Step ===================================== - def step_hyperbolic(self,solution): - r""" - Take a step on the homogeneous hyperbolic system using the Clawpack - algorithm. - - Clawpack is based on the Lax-Wendroff method, combined with Riemann - solvers and TVD limiters applied to waves. - """ - if(self.kernel_language == 'Fortran'): - state = solution.states[0] - grid = state.grid - dx,dy,dz = grid.delta - mx,my,mz = grid.num_cells - maxm = max(mx,my,mz) - - self._apply_bcs(state) - qnew = self.qbc - qold = qnew.copy('F') - - rpn3 = self.rp.rpn3._cpointer - - if (self.dimensional_split) or (self.transverse_waves==0): - rpt3 = rpn3 # dummy value; it won't be called - rptt3 = rpn3 # dummy value; it won't be called - else: - rpt3 = self.rp.rpt3._cpointer - rptt3 = self.rp.rptt3._cpointer - - if self.dimensional_split: - #Right now only Godunov-dimensional-splitting is implemented. - #Strang-dimensional-splitting could be added following dimsp3.f in Clawpack. - - q, cfl_x = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ - qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,1,self.fwave,rpn3,rpt3,rptt3) - - q, cfl_y = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ - q,q,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,2,self.fwave,rpn3,rpt3,rptt3) - - q, cfl_z = self.fmod.step3ds(maxm,self.num_ghost,mx,my,mz, \ - q,q,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,3,self.fwave,rpn3,rpt3,rptt3) - - cfl = max(cfl_x,cfl_y,cfl_z) - - else: - - q, cfl = self.fmod.step3(maxm,self.num_ghost,mx,my,mz, \ - qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,self._mthlim,\ - self.aux1,self.aux2,self.aux3,self.work,self.fwave,rpn3,rpt3,rptt3) - - self.cfl.update_global_max(cfl) - state.set_q_from_qbc(self.num_ghost,self.qbc) - if state.num_aux > 0: - state.set_aux_from_auxbc(self.num_ghost,self.auxbc) - - else: - raise NotImplementedError("No python implementation for step_hyperbolic in 3D.")