Merge 4a25c38 into b375bcf

examples/advection_1d/advection_1d_nonunif.py

@@ -0,0 +1,121 @@
+#!/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
+
+    # sharpclaw does not accomodate nonuniform grids
+    if solver_type=='classic':
+        solver = pyclaw.ClawSolver1D(riemann_solver)
+    else: raise Exception('Unrecognized value of solver_type.')
+
+    solver.kernel_language = kernel_language
+    solver.order = 1 
+    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)

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()

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]
 

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: ]