From ae421626555281aae6c1e4fcd695448e7c60a503 Mon Sep 17 00:00:00 2001 From: cr2940 Date: Tue, 22 Oct 2019 18:18:07 -0400 Subject: [PATCH 1/2] added nonuniform examples and fixed aux array bug in classic and sharpclaw --- examples/advection_1d/advection_1d_nonunif.py | 127 ++++++++++++++++++ .../advection_1d/test_advection_nonunif.py | 83 ++++++++++++ src/pyclaw/classic/solver.py | 2 +- src/pyclaw/sharpclaw/solver.py | 6 +- 4 files changed, 214 insertions(+), 4 deletions(-) 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..0188aa738 --- /dev/null +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -0,0 +1,127 @@ +#!/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. +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 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 [-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] + """ + neg = -1*(xc < 0) + (xc > 0) + xp = xc**2 + xp = neg*xp + return xp + + +def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='classic', weno_order=5, + time_integrator='SSP104', outdir='./_output'): + + if use_petsc: + import clawpack.petclaw as pyclaw + from clawpack import pyclaw + import clawpack.pyclaw.geometry + + 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 # in order to make test file pass without complaint + solver.limiters = pyclaw.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 = 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)) # 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.0 + 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 = 0.5 # one cycle + claw.outdir = outdir + claw.num_output_times = 10 + claw.nstepout = 1 + if outdir is None: + claw.output_format = None + + claw.setplot = setplot + + 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 = [-0.25,0.25] + 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) diff --git a/examples/advection_1d/test_advection_nonunif.py b/examples/advection_1d/test_advection_nonunif.py new file mode 100644 index 000000000..7b6936233 --- /dev/null +++ b/examples/advection_1d/test_advection_nonunif.py @@ -0,0 +1,83 @@ +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 + + 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 + + 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 + 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))) + 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_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.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.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(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, sharp_tests_lmm2, weno_tests): + yield test + + +if __name__=="__main__": + import nose + nose.main() 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] diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 02b90f5bd..1c5401765 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -735,7 +735,8 @@ def dq_hyperbolic(self,state): import numpy as np self._apply_bcs(state) - q = self.qbc + q = self.qbc + aux = self.auxbc grid = state.grid mx = grid.num_cells[0] @@ -758,11 +759,10 @@ def dq_hyperbolic(self,state): # 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: ] From 4a25c38d34e3cfb81a87db89e41df5ee65ffc6b1 Mon Sep 17 00:00:00 2001 From: cr2940 <36740525+cr2940@users.noreply.github.com> Date: Wed, 23 Oct 2019 12:26:34 -0400 Subject: [PATCH 2/2] remove sharpclaw solver option --- examples/advection_1d/advection_1d_nonunif.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 0188aa738..3df68c528 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/examples/advection_1d/advection_1d_nonunif.py @@ -45,19 +45,13 @@ def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='class elif kernel_language == 'Python': riemann_solver = riemann.advection_1D_py.advection_1D + # sharpclaw does not accomodate nonuniform grids 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 # in order to make test file pass without complaint + solver.order = 1 solver.limiters = pyclaw.tvd.minmod solver.num_eqn=1 solver.num_waves=1