Polynomial reconstruction in SharpClaw

examples/acoustics_1d_homogeneous/acoustics_1d.py

@@ -24,7 +24,7 @@
 
 
 def setup(use_petsc=False, kernel_language='Fortran', solver_type='classic',
-          outdir='./_output', ptwise=False, weno_order=5,
+          outdir='./_output', ptwise=False, reconstruction_order=5, lim_type=2,
           time_integrator='SSP104', disable_output=False, output_style=1):
 
     if use_petsc:
@@ -46,8 +46,9 @@ def setup(use_petsc=False, kernel_language='Fortran', solver_type='classic',
         solver.limiters = pyclaw.limiters.tvd.MC
     elif solver_type == 'sharpclaw':
         solver = pyclaw.SharpClawSolver1D(riemann_solver)
-        solver.weno_order = weno_order
+        solver.reconstruction_order = reconstruction_order
         solver.time_integrator = time_integrator
+        solver.lim_type = lim_type
         if time_integrator == 'SSPLMMk3':
             solver.lmm_steps = 4
     else:

examples/acoustics_1d_homogeneous/test_acoustics.py

@@ -59,7 +59,7 @@ def acoustics_verify(claw):
     weno_tests = gen_variants(acoustics_1d.setup, verify_expected(0.000153),
                               kernel_languages=('Fortran',),
                               solver_type='sharpclaw', time_integrator='SSP104',
-                              weno_order=17, disable_output=True)
+                              reconstruction_order=17, disable_output=True)
 
     from itertools import chain
     for test in chain(classic_tests, time_step_test, ptwise_tests,

examples/advection_1d/advection_1d.py

@@ -20,7 +20,7 @@
 import numpy as np
 from clawpack import riemann
 
-def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, 
+def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classic', reconstruction_order=5, 
           time_integrator='SSP104', outdir='./_output'):
 
     if use_petsc:
@@ -37,7 +37,7 @@ def setup(nx=100, kernel_language='Python', use_petsc=False, solver_type='classi
         solver = pyclaw.ClawSolver1D(riemann_solver)
     elif solver_type=='sharpclaw':
         solver = pyclaw.SharpClawSolver1D(riemann_solver)
-        solver.weno_order = weno_order
+        solver.reconstruction_order = reconstruction_order
         solver.time_integrator = time_integrator
         if time_integrator == 'SSPLMMk3':
             solver.lmm_steps = 5

examples/advection_1d/advection_1d_nonunif.py

@@ -32,7 +32,7 @@ def mapc2p_nonunif(xc):
     return xp
 
 
-def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='classic', weno_order=5,
+def setup(nx=100, kernel_language='Fortran', use_petsc=False, solver_type='classic', reconstruction_order=5,
           time_integrator='SSP104', outdir='./_output'):
 
     if use_petsc:

examples/advection_1d/test_advection.py

@@ -39,7 +39,7 @@ def advection_verify(claw):
 
     weno_tests = gen_variants(advection_1d.setup, verify_expected(7.489618e-06),
                                  kernel_languages=('Fortran',), solver_type='sharpclaw', 
-                                 time_integrator='SSP104', weno_order=17,
+                                 time_integrator='SSP104', reconstruction_order=17,
                                  outdir=None)
 
     from itertools import chain

examples/advection_1d_variable/variable_coefficient_advection.py

@@ -69,7 +69,7 @@ def setup(use_petsc=False,solver_type='classic',kernel_language='Python',outdir=
             solver = pyclaw.SharpClawSolver1D(riemann.advection_color_1D)
         elif kernel_language=='Python': 
             solver = pyclaw.SharpClawSolver1D(riemann.vc_advection_1D_py.vc_advection_1D)
-        solver.weno_order=weno_order
+        solver.reconstruction_order=reconstruction_order
     else: raise Exception('Unrecognized value of solver_type.')
 
     solver.kernel_language = kernel_language

examples/cubic_1d/cubic.py

@@ -17,7 +17,7 @@
 import numpy as np
 from clawpack import riemann
 
-def setup(use_petsc=0, outdir='./_output', solver_type='classic', weno_order=5, N=1000):
+def setup(use_petsc=0, outdir='./_output', solver_type='classic', reconstruction_order=5, lim_type=2, N=1000):
 
     if use_petsc:
         import clawpack.petclaw as pyclaw
@@ -28,7 +28,8 @@ def setup(use_petsc=0, outdir='./_output', solver_type='classic', weno_order=5,
 
     if solver_type=='sharpclaw':
         solver = pyclaw.SharpClawSolver1D(riemann_solver)
-        solver.weno_order = weno_order
+        solver.reconstruction_order = reconstruction_order
+        solver.lim_type = lim_type
     else:
         solver = pyclaw.ClawSolver1D(riemann_solver)
         solver.limiters = pyclaw.limiters.tvd.vanleer

examples/euler_1d/setup.py

@@ -22,6 +22,7 @@ def configuration(parent_package='',top_path=None):
     sharpclaw_srcs = [pjoin(sharpclaw_dir, src) for src in ['ClawParams.f90',
                                                             'workspace.f90',
                                                             'weno.f90',
+                                                            'poly.f90',
                                                             'reconstruct.f90','flux1.f90']]
 
     config.add_extension('sharpclaw1',

examples/euler_2d/shock_bubble_interaction.py

@@ -186,7 +186,7 @@ def setup(use_petsc=False,solver_type='classic', outdir='_output', kernel_langua
     if solver_type=='sharpclaw':
         solver = pyclaw.SharpClawSolver2D(riemann.euler_5wave_2D)
         solver.dq_src = dq_Euler_radial
-        solver.weno_order = 5
+        solver.reconstruction_order = 5
         solver.lim_type   = 2
     else:
         solver = pyclaw.ClawSolver2D(riemann.euler_5wave_2D)

fvmbook/chap6/limiters.ipynb

@@ -55,7 +55,7 @@
     "    elif scheme == 'first-order':\n",
     "        solver.order = 1\n",
     "    elif 'weno' in scheme:\n",
-    "        solver.weno_order = int(scheme[4:]) #weno5, weno7, ...\n",
+    "        solver.reconstruction_order = int(scheme[4:]) #weno5, weno7, ...\n",
     "    else:\n",
     "        raise Exception('Unrecognized limiter')\n",
     "\n",

src/pyclaw/sharpclaw/ClawParams.f90

@@ -9,7 +9,7 @@ module ClawParams
     integer :: num_dim, num_waves, index_capa
 
     ! Method-related parameters:
-    integer :: char_decomp, lim_type, weno_order
+    integer :: char_decomp, lim_type, reconstruction_order
     integer, allocatable :: mthlim(:)
     logical :: fwave, tfluct_solver
 

src/pyclaw/sharpclaw/flux1.f90

@@ -115,9 +115,11 @@ subroutine flux1(q1d,dq1d,aux,dt,cfl,t,ixyz,num_aux,num_eqn,mx,num_ghost,maxnx,r
         end select
         case(3)
             call weno5(q1d,ql,qr,num_eqn,maxnx,num_ghost)
+        case(4)
+            call poly_comp(q1d,ql,qr,num_eqn,maxnx,num_ghost)
         case default
             write(*,*) 'ERROR: Unrecognized limiter type option'
-            write(*,*) 'You should set 1<=lim_type<=3'
+            write(*,*) 'You should set 1<=lim_type<=4'
             stop
       end select
 

src/pyclaw/sharpclaw/poly.f90

@@ -0,0 +1,95 @@
+module poly
+contains
+
+    ! ===================================================================
+    subroutine poly4(q,ql,qr,num_eqn,maxnx,num_ghost)
+    ! ===================================================================
+
+        implicit none
+
+        integer,          intent(in) :: num_eqn, maxnx, num_ghost
+        double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost)
+        double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost)
+
+        integer :: i,n
+
+        do n = 1,num_eqn
+            do i = num_ghost-1,maxnx+num_ghost+1
+                ql(n,i) = (-5.d0*q(n,-2+i)+60.d0*q(n,-1+i)+90.d0*q(n,i)-20.d0*q(n,1+i)+3.d0*q(n,2+i))/128.d0
+                qr(n,i) = (3.d0*q(n,-2+i)-20.d0*q(n,-1+i)+90.d0*q(n,i)+60.d0*q(n,1+i)-5.d0*q(n,2+i))/128.d0
+            end do
+        end do
+
+    end subroutine poly4
+
+    ! ===================================================================
+    subroutine poly6(q,ql,qr,num_eqn,maxnx,num_ghost)
+    ! ===================================================================
+
+        implicit none
+
+        integer,          intent(in) :: num_eqn, maxnx, num_ghost
+        double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost)
+        double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost)
+
+        integer :: i,n
+
+        do n = 1,num_eqn
+            do i = num_ghost-1,maxnx+num_ghost+1
+                ql(n,i) = (7.d0*q(n,-3+i) - 70.d0*q(n,-2+i) + 525.d0*q(n,-1+i) + 700.d0*q(n,i) &
+                           - 175.d0*q(n,1+i) + 42.d0*q(n,2+i) - 5.d0*q(n,3+i))/1024.d0
+                qr(n,i) = (42.d0*q(n,-2+i) - 175.d0*q(n,-1+i) + 700.d0*q(n,i) &
+                           + 525.d0*q(n,1+i) - 70.d0*q(n,2+i) + 7.d0*q(n,3+i) - 5.d0*q(n,-3+i))/1024.d0
+            end do
+        end do
+
+    end subroutine poly6
+
+    ! ===================================================================
+    subroutine poly8(q,ql,qr,num_eqn,maxnx,num_ghost)
+    ! ===================================================================
+
+        implicit none
+
+        integer,          intent(in) :: num_eqn, maxnx, num_ghost
+        double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost)
+        double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost)
+
+        integer :: i,n
+
+        do n = 1,num_eqn
+            do i = num_ghost-1,maxnx+num_ghost+1
+                ql(n,i) = (-45.d0*q(n,-4+i) + 504.d0*q(n,-3+i) - 2940.d0*q(n,-2+i) + 17640.d0*q(n,-1+i) &
+                           + 22050.d0*q(n,i) - 5880.d0*q(n,1+i) + 1764.d0*q(n,2+i) - 360.d0*q(n,3+i) + 35.d0*q(n,4+i))/32768.d0
+                qr(n,i) = (35.d0*q(n,-4+i) - 360.d0*q(n,-3+i) + 1764.d0*q(n,-2+i) - 5880.d0*q(n,-1+i) + 22050.d0*q(n,i) &
+                           + 17640.d0*q(n,1+i) - 2940.d0*q(n,2+i) + 504.d0*q(n,3+i) - 45.d0*q(n,4+i))/32768.d0
+            end do
+        end do
+
+    end subroutine poly8
+
+    ! ===================================================================
+    subroutine poly10(q,ql,qr,num_eqn,maxnx,num_ghost)
+    ! ===================================================================
+
+        implicit none
+
+        integer,          intent(in) :: num_eqn, maxnx, num_ghost
+        double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost)
+        double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost)
+
+        integer :: i,n
+
+        do n = 1,num_eqn
+            do i = num_ghost-1,maxnx+num_ghost+1
+                ql(n,i) = (77.d0*q(n,-5+i) - 990.d0*q(n,-4+i) + 6237.d0*q(n,-3+i) - 27720.d0*q(n,-2+i) &
+                           + 145530.d0*q(n,-1+i) + 174636.d0*q(n,i) - 48510.d0*q(n,1+i) + 16632.d0*q(n,2+i) - 4455.d0*q(n,3+i) &
+                           + 770.d0*q(n,4+i) - 63.d0*q(n,5+i))/262144.d0
+                qr(n,i) = (-63.d0*q(n,-5+i) + 770.d0*q(n,-4+i) - 4455.d0*q(n,-3+i) + 16632.d0*q(n,-2+i) & 
+                           - 48510.d0*q(n,-1+i) + 174636.d0*q(n,i) + 145530.d0*q(n,1+i) - 27720.d0*q(n,2+i) + 6237.d0*q(n,3+i) &
+                           - 990.d0*q(n,4+i) + 77.d0*q(n,5+i))/262144.d0
+            end do
+        end do
+    end subroutine poly10
+
+end module poly

src/pyclaw/sharpclaw/reconstruct.f90

@@ -83,8 +83,12 @@ subroutine dealloc_recon_workspace(lim_type,char_decomp)
                     deallocate(hh)
                     deallocate(uh)
             end select
-            recon_alloc = .False.
+            case(3)
+                deallocate(uu)
+                deallocate(dq1m)
         end select
+        recon_alloc = .False.
+
     end subroutine dealloc_recon_workspace
 
 ! ===================================================================
@@ -97,14 +101,14 @@ subroutine weno_comp(q,ql,qr,num_eqn,maxnx,num_ghost)
     !   It does no characteristic decomposition
 
         use weno
-        use clawparams, only: weno_order
+        use clawparams, only: reconstruction_order
         implicit none
 
         integer,          intent(in) :: num_eqn, maxnx, num_ghost
         double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost)
         double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost)
 
-        select case(weno_order)
+        select case(reconstruction_order)
         case (5)
            call weno5(q,ql,qr,num_eqn,maxnx,num_ghost)
         case (7)
@@ -120,7 +124,7 @@ subroutine weno_comp(q,ql,qr,num_eqn,maxnx,num_ghost)
         case (17)
            call weno17(q,ql,qr,num_eqn,maxnx,num_ghost)           
         case default
-           print *, 'ERROR: weno_order must be an odd number between 5 and 17 (inclusive).'
+           print *, 'ERROR: reconstruction_order in WENO must be an odd number between 5 and 17 (inclusive).'
            stop
         end select
 
@@ -821,4 +825,35 @@ subroutine tvd2_wave(q,ql,qr,wave,s,mthlim,num_eqn,num_ghost)
       return
       end subroutine tvd2_wave
 
+    ! ===================================================================
+    subroutine poly_comp(q,ql,qr,num_eqn,maxnx,num_ghost)
+    ! ===================================================================
+
+        use poly
+        use clawparams, only: reconstruction_order
+        implicit none
+
+        integer,          intent(in) :: num_eqn, maxnx, num_ghost
+        double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost)
+        double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost)
+
+        select case(reconstruction_order)
+        case (4)
+            call poly4(q,ql,qr,num_eqn,maxnx,num_ghost)
+        case (6)
+            call poly6(q,ql,qr,num_eqn,maxnx,num_ghost)           
+        case (8)
+            print *, 'Warning: 8th reconstruction order is not fully tested'
+            call poly8(q,ql,qr,num_eqn,maxnx,num_ghost)
+        case(10)
+            print *, 'Warning: 10th reconstruction order is not fully tested'
+            call poly10(q,ql,qr,num_eqn,maxnx,num_ghost)          
+        case default
+           print *, 'ERROR: reconstruction_order for polynomial must be an odd number between 4 and 10 (inclusive)'
+           stop
+        end select
+
+    return
+    end subroutine poly_comp
+
 end module reconstruct

src/pyclaw/sharpclaw/setup.py

@@ -6,16 +6,16 @@ def configuration(parent_package='',top_path=None):
     config = Configuration('sharpclaw', parent_package, top_path)
 
     config.add_extension('sharpclaw1',
-                         ['ClawParams.f90','weno.f90','reconstruct.f90',
+                         ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90',
                           'evec.f90','workspace.f90','flux1.f90'])
 
     config.add_extension('sharpclaw2',
-                         ['ClawParams.f90','weno.f90','reconstruct.f90',
+                         ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90',
                           'evec.f90','workspace.f90','flux2.f90',
                           'flux1.f90'])
 
     config.add_extension('sharpclaw3',
-                         ['ClawParams.f90','weno.f90','reconstruct.f90',
+                         ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90',
                           'evec.f90','workspace.f90','flux3.f90',
                           'flux1.f90'])
 

src/pyclaw/sharpclaw/solver.py

@@ -42,13 +42,14 @@ class SharpClawSolver(Solver):
             - 0: No limiting.
             - 1: TVD reconstruction.
             - 2: WENO reconstruction.
+            - 3: Polynomial reconstruction.
 
         ``Default = 2``
 
-    .. attribute:: weno_order
-
-        Order of the WENO reconstruction. From 1st to 17th order (PyWENO)
+    .. attribute:: reconstruction_order
 
+        Order of the reconstruction. From 1st to 17th order (PyWENO), and from
+        4 to 10 (even) in Poly
         ``Default = 5``
 
     .. attribute:: time_integrator
@@ -161,7 +162,7 @@ def __init__(self,riemann_solver=None,claw_package=None):
         """
         self.limiters = [1]
         self.lim_type = 2
-        self.weno_order = 5
+        self.reconstruction_order = 5
         self.time_integrator = 'SSP104'
         self.char_decomp = 0
         self.tfluct_solver = False
@@ -206,9 +207,12 @@ def setup(self,solution):
         Allocate RK stage arrays or previous step solutions and fortran routine work arrays.
         """
         if self.lim_type == 2:
-            self.num_ghost = (self.weno_order+1)//2
+            self.num_ghost = int((self.reconstruction_order+1)/2)
+
+        if self.lim_type == 3:
+            self.num_ghost = int(1 + self.reconstruction_order/2)
 
-        if self.lim_type == 2 and self.weno_order != 5 and self.kernel_language == 'Python':
+        if self.lim_type == 2 and self.reconstruction_order != 5 and self.kernel_language == 'Python':
             raise Exception('Only 5th-order WENO reconstruction is implemented in Python kernels. \
                              Use Fortran for higher-order WENO.')
         # This is a hack to deal with the fact that petsc4py
@@ -565,7 +569,7 @@ def _set_fortran_parameters(self,state,clawparams,workspace,reconstruct):
         grid = state.grid
         clawparams.num_dim       = grid.num_dim
         clawparams.lim_type      = self.lim_type
-        clawparams.weno_order    = self.weno_order
+        clawparams.reconstruction_order    = self.reconstruction_order
         clawparams.char_decomp   = self.char_decomp
         clawparams.tfluct_solver = self.tfluct_solver
         clawparams.fwave         = self.fwave