From 1d44160d58a66b10dffc0d8b648f82739edd35c1 Mon Sep 17 00:00:00 2001 From: sanromd Date: Tue, 2 Dec 2014 16:43:25 +0300 Subject: [PATCH 01/14] polynomial interpolation routines added --- src/pyclaw/sharpclaw/poly.f90 | 95 +++++++++++++++++++++++++++++++++++ 1 file changed, 95 insertions(+) create mode 100644 src/pyclaw/sharpclaw/poly.f90 diff --git a/src/pyclaw/sharpclaw/poly.f90 b/src/pyclaw/sharpclaw/poly.f90 new file mode 100644 index 000000000..cb095a7a3 --- /dev/null +++ b/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 \ No newline at end of file From 62ecd182e4688495d86f1ecd25d38ebb15fcf22c Mon Sep 17 00:00:00 2001 From: sanromd Date: Tue, 25 Jul 2017 16:10:04 +0300 Subject: [PATCH 02/14] solved merge conflicts by changing weno_order to interpolation_order --- examples/acoustics_1d_homogeneous/acoustics_1d.py | 4 ++-- examples/acoustics_1d_homogeneous/test_acoustics.py | 2 +- examples/advection_1d/advection_1d.py | 4 ++-- examples/advection_1d/test_advection.py | 2 +- .../variable_coefficient_advection.py | 2 +- examples/euler_2d/shock_bubble_interaction.py | 2 +- src/pyclaw/sharpclaw/ClawParams.f90 | 2 +- src/pyclaw/sharpclaw/reconstruct.f90 | 6 +++--- src/pyclaw/sharpclaw/solver.py | 10 +++++----- 9 files changed, 17 insertions(+), 17 deletions(-) diff --git a/examples/acoustics_1d_homogeneous/acoustics_1d.py b/examples/acoustics_1d_homogeneous/acoustics_1d.py index 72f1f0e39..2464ffa90 100644 --- a/examples/acoustics_1d_homogeneous/acoustics_1d.py +++ b/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, interpolation_order=5, time_integrator='SSP104', disable_output=False, output_style=1): if use_petsc: @@ -46,7 +46,7 @@ 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.interpolation_order = interpolation_order solver.time_integrator = time_integrator if time_integrator == 'SSPLMMk3': solver.lmm_steps = 4 diff --git a/examples/acoustics_1d_homogeneous/test_acoustics.py b/examples/acoustics_1d_homogeneous/test_acoustics.py index 83682a054..0c5e3bdef 100644 --- a/examples/acoustics_1d_homogeneous/test_acoustics.py +++ b/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) + interpolation_order=17, disable_output=True) from itertools import chain for test in chain(classic_tests, time_step_test, ptwise_tests, diff --git a/examples/advection_1d/advection_1d.py b/examples/advection_1d/advection_1d.py index 533505165..c6d624b33 100755 --- a/examples/advection_1d/advection_1d.py +++ b/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', interpolation_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.interpolation_order = interpolation_order solver.time_integrator = time_integrator if time_integrator == 'SSPLMMk3': solver.lmm_steps = 5 diff --git a/examples/advection_1d/test_advection.py b/examples/advection_1d/test_advection.py index edd54bb8a..8b5fe3d86 100644 --- a/examples/advection_1d/test_advection.py +++ b/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', interpolation_order=17, outdir=None) from itertools import chain diff --git a/examples/advection_1d_variable/variable_coefficient_advection.py b/examples/advection_1d_variable/variable_coefficient_advection.py index 4dd01af25..89f80dc50 100755 --- a/examples/advection_1d_variable/variable_coefficient_advection.py +++ b/examples/advection_1d_variable/variable_coefficient_advection.py @@ -68,7 +68,7 @@ def setup(use_petsc=False,solver_type='classic',kernel_language='Python',outdir= solver = pyclaw.SharpClawSolver1D(riemann.vc_advection_1D) elif kernel_language=='Python': solver = pyclaw.SharpClawSolver1D(riemann.vc_advection_1D_py.vc_advection_1D) - solver.weno_order=weno_order + solver.interpolation_order=interpolation_order else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language diff --git a/examples/euler_2d/shock_bubble_interaction.py b/examples/euler_2d/shock_bubble_interaction.py index 6d847bb3f..d97b72e8e 100755 --- a/examples/euler_2d/shock_bubble_interaction.py +++ b/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.interpolation_order = 5 solver.lim_type = 2 else: solver = pyclaw.ClawSolver2D(riemann.euler_5wave_2D) diff --git a/src/pyclaw/sharpclaw/ClawParams.f90 b/src/pyclaw/sharpclaw/ClawParams.f90 index b1e28357e..81d8a6764 100644 --- a/src/pyclaw/sharpclaw/ClawParams.f90 +++ b/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, interpolation_order integer, allocatable :: mthlim(:) logical :: fwave, tfluct_solver diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 4db69532e..6502ec8dc 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -97,14 +97,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: interpolation_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(interpolation_order) case (5) call weno5(q,ql,qr,num_eqn,maxnx,num_ghost) case (7) @@ -120,7 +120,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: interpolation_order must be an odd number between 5 and 17 (inclusive).' stop end select diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 1c5401765..4e96cfe79 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -45,7 +45,7 @@ class SharpClawSolver(Solver): ``Default = 2`` - .. attribute:: weno_order + .. attribute:: interpolation_order Order of the WENO reconstruction. From 1st to 17th order (PyWENO) @@ -161,7 +161,7 @@ def __init__(self,riemann_solver=None,claw_package=None): """ self.limiters = [1] self.lim_type = 2 - self.weno_order = 5 + self.interpolation_order = 5 self.time_integrator = 'SSP104' self.char_decomp = 0 self.tfluct_solver = False @@ -206,9 +206,9 @@ 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 = (self.interpolation_order+1)/2 - if self.lim_type == 2 and self.weno_order != 5 and self.kernel_language == 'Python': + if self.lim_type == 2 and self.interpolation_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 @@ -563,7 +563,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.interpolation_order = self.interpolation_order clawparams.char_decomp = self.char_decomp clawparams.tfluct_solver = self.tfluct_solver clawparams.fwave = self.fwave From f9edcb104fd0dcff159fe447a3d22980a40d1802 Mon Sep 17 00:00:00 2001 From: sanromd Date: Wed, 3 Dec 2014 14:59:34 +0300 Subject: [PATCH 03/14] added calls to module poly in resonsctruct and set setup to compile poly --- src/pyclaw/sharpclaw/reconstruct.f90 | 29 ++++++++++++++++++++++++++++ src/pyclaw/sharpclaw/setup.py | 6 +++--- 2 files changed, 32 insertions(+), 3 deletions(-) diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 6502ec8dc..91378cf74 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -821,4 +821,33 @@ 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: interpolation_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(interpolation_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) + call poly8(q,ql,qr,num_eqn,maxnx,num_ghost) + case(10) + call poly10(q,ql,qr,num_eqn,maxnx,num_ghost) + case default + print *, 'ERROR: poly_order must be an even number between 4 and 8.' + stop + end select + + return + end subroutine poly_comp + end module reconstruct diff --git a/src/pyclaw/sharpclaw/setup.py b/src/pyclaw/sharpclaw/setup.py index 38c076780..d96369e3a 100644 --- a/src/pyclaw/sharpclaw/setup.py +++ b/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']) From ffe513fff9ec309f5b27492c1d97b3106d604afc Mon Sep 17 00:00:00 2001 From: sanromd Date: Wed, 3 Dec 2014 17:25:51 +0300 Subject: [PATCH 04/14] added WARNING and ERROR messages to reconsctruct when reconstrcution order for Poly is >6 --- src/pyclaw/sharpclaw/reconstruct.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 91378cf74..5e1d5ae20 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -120,7 +120,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: interpolation_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 @@ -839,11 +839,13 @@ subroutine poly_comp(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: poly_order must be an even number between 4 and 8.' + print *, 'ERROR: reconstruction_order for polynomial must be either 4 or 6' stop end select From 9a648ee661dbe6cdcfe0a275d3aac4c90dc6743e Mon Sep 17 00:00:00 2001 From: sanromd Date: Wed, 3 Dec 2014 17:26:14 +0300 Subject: [PATCH 05/14] added poly to custom sharpclaw build in euler_1d --- examples/euler_1d/setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/euler_1d/setup.py b/examples/euler_1d/setup.py index 0bc16b985..6cdbbf5dd 100644 --- a/examples/euler_1d/setup.py +++ b/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', From 1dbfabe909fcabad46e350d393808e7b5f3ace79 Mon Sep 17 00:00:00 2001 From: sanromd Date: Tue, 25 Jul 2017 16:12:29 +0300 Subject: [PATCH 06/14] solved merge conflict by changin weno_order to interpolation_order --- examples/acoustics_1d_homogeneous/acoustics_1d.py | 4 ++-- examples/acoustics_1d_homogeneous/test_acoustics.py | 2 +- examples/advection_1d/advection_1d.py | 4 ++-- examples/advection_1d/test_advection.py | 2 +- .../variable_coefficient_advection.py | 2 +- examples/euler_2d/shock_bubble_interaction.py | 2 +- src/pyclaw/sharpclaw/ClawParams.f90 | 2 +- src/pyclaw/sharpclaw/reconstruct.f90 | 8 ++++---- src/pyclaw/sharpclaw/solver.py | 10 +++++----- 9 files changed, 18 insertions(+), 18 deletions(-) diff --git a/examples/acoustics_1d_homogeneous/acoustics_1d.py b/examples/acoustics_1d_homogeneous/acoustics_1d.py index 2464ffa90..b06187783 100644 --- a/examples/acoustics_1d_homogeneous/acoustics_1d.py +++ b/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, interpolation_order=5, + outdir='./_output', ptwise=False, reconstruction_order=5, time_integrator='SSP104', disable_output=False, output_style=1): if use_petsc: @@ -46,7 +46,7 @@ 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.interpolation_order = interpolation_order + solver.reconstruction_order = reconstruction_order solver.time_integrator = time_integrator if time_integrator == 'SSPLMMk3': solver.lmm_steps = 4 diff --git a/examples/acoustics_1d_homogeneous/test_acoustics.py b/examples/acoustics_1d_homogeneous/test_acoustics.py index 0c5e3bdef..61e53245c 100644 --- a/examples/acoustics_1d_homogeneous/test_acoustics.py +++ b/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', - interpolation_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, diff --git a/examples/advection_1d/advection_1d.py b/examples/advection_1d/advection_1d.py index c6d624b33..2e68c1306 100755 --- a/examples/advection_1d/advection_1d.py +++ b/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', interpolation_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.interpolation_order = interpolation_order + solver.reconstruction_order = reconstruction_order solver.time_integrator = time_integrator if time_integrator == 'SSPLMMk3': solver.lmm_steps = 5 diff --git a/examples/advection_1d/test_advection.py b/examples/advection_1d/test_advection.py index 8b5fe3d86..53c35f49d 100644 --- a/examples/advection_1d/test_advection.py +++ b/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', interpolation_order=17, + time_integrator='SSP104', reconstruction_order=17, outdir=None) from itertools import chain diff --git a/examples/advection_1d_variable/variable_coefficient_advection.py b/examples/advection_1d_variable/variable_coefficient_advection.py index 89f80dc50..347ba76df 100755 --- a/examples/advection_1d_variable/variable_coefficient_advection.py +++ b/examples/advection_1d_variable/variable_coefficient_advection.py @@ -68,7 +68,7 @@ def setup(use_petsc=False,solver_type='classic',kernel_language='Python',outdir= solver = pyclaw.SharpClawSolver1D(riemann.vc_advection_1D) elif kernel_language=='Python': solver = pyclaw.SharpClawSolver1D(riemann.vc_advection_1D_py.vc_advection_1D) - solver.interpolation_order=interpolation_order + solver.reconstruction_order=reconstruction_order else: raise Exception('Unrecognized value of solver_type.') solver.kernel_language = kernel_language diff --git a/examples/euler_2d/shock_bubble_interaction.py b/examples/euler_2d/shock_bubble_interaction.py index d97b72e8e..934cf2dc9 100755 --- a/examples/euler_2d/shock_bubble_interaction.py +++ b/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.interpolation_order = 5 + solver.reconstruction_order = 5 solver.lim_type = 2 else: solver = pyclaw.ClawSolver2D(riemann.euler_5wave_2D) diff --git a/src/pyclaw/sharpclaw/ClawParams.f90 b/src/pyclaw/sharpclaw/ClawParams.f90 index 81d8a6764..75f671711 100644 --- a/src/pyclaw/sharpclaw/ClawParams.f90 +++ b/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, interpolation_order + integer :: char_decomp, lim_type, reconstruction_order integer, allocatable :: mthlim(:) logical :: fwave, tfluct_solver diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 5e1d5ae20..e9cf7df80 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -97,14 +97,14 @@ subroutine weno_comp(q,ql,qr,num_eqn,maxnx,num_ghost) ! It does no characteristic decomposition use weno - use clawparams, only: interpolation_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(interpolation_order) + select case(reconstruction_order) case (5) call weno5(q,ql,qr,num_eqn,maxnx,num_ghost) case (7) @@ -826,14 +826,14 @@ subroutine poly_comp(q,ql,qr,num_eqn,maxnx,num_ghost) ! =================================================================== use poly - use clawparams, only: interpolation_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(interpolation_order) + select case(reconstruction_order) case (4) call poly4(q,ql,qr,num_eqn,maxnx,num_ghost) case (6) diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 4e96cfe79..c71a34919 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -45,7 +45,7 @@ class SharpClawSolver(Solver): ``Default = 2`` - .. attribute:: interpolation_order + .. attribute:: reconstruction_order Order of the WENO reconstruction. From 1st to 17th order (PyWENO) @@ -161,7 +161,7 @@ def __init__(self,riemann_solver=None,claw_package=None): """ self.limiters = [1] self.lim_type = 2 - self.interpolation_order = 5 + self.reconstruction_order = 5 self.time_integrator = 'SSP104' self.char_decomp = 0 self.tfluct_solver = False @@ -206,9 +206,9 @@ 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.interpolation_order+1)/2 + self.num_ghost = (self.reconstruction_order+1)/2 - if self.lim_type == 2 and self.interpolation_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 @@ -563,7 +563,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.interpolation_order = self.interpolation_order + clawparams.reconstruction_order = self.reconstruction_order clawparams.char_decomp = self.char_decomp clawparams.tfluct_solver = self.tfluct_solver clawparams.fwave = self.fwave From 5a53304020fe87cc542a46d9c13e7d43cb68f96c Mon Sep 17 00:00:00 2001 From: sanromd Date: Sun, 28 Dec 2014 14:13:36 +0300 Subject: [PATCH 07/14] added allocated, deallocate and call for polyinterpolation and WENO5 --- src/pyclaw/sharpclaw/flux1.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pyclaw/sharpclaw/flux1.f90 b/src/pyclaw/sharpclaw/flux1.f90 index e631adaa5..7e9ccce36 100644 --- a/src/pyclaw/sharpclaw/flux1.f90 +++ b/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 From 943c2f2c0bab850668d77651260ef92b2d59947f Mon Sep 17 00:00:00 2001 From: sanromd Date: Sun, 28 Dec 2014 14:13:56 +0300 Subject: [PATCH 08/14] allocate in reconstruct --- src/pyclaw/sharpclaw/reconstruct.f90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index e9cf7df80..781b866fb 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -48,6 +48,9 @@ subroutine alloc_recon_workspace(maxnx,num_ghost,num_eqn,num_waves,lim_type,char case(3) allocate(uu(2,maxnx+2*num_ghost)) allocate( dq1m(maxnx+2*num_ghost)) + case(4) + allocate(uu(2,maxnx+2*num_ghost)) + allocate( dq1m(maxnx+2*num_ghost)) end select recon_alloc = .True. @@ -83,8 +86,15 @@ subroutine dealloc_recon_workspace(lim_type,char_decomp) deallocate(hh) deallocate(uh) end select - recon_alloc = .False. + case(3) + deallocate(uu) + deallocate(dq1m) + case(4) + deallocate(uu) + deallocate(dq1m) end select + recon_alloc = .False. + end subroutine dealloc_recon_workspace ! =================================================================== From f7f409f831f651556fd836ff8345a8215ee033ff Mon Sep 17 00:00:00 2001 From: sanromd Date: Sun, 28 Dec 2014 14:14:31 +0300 Subject: [PATCH 09/14] added documentation to solver for lim_type=4 (polynomial) --- src/pyclaw/sharpclaw/solver.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index c71a34919..cfcb4e278 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/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:: reconstruction_order - Order of the WENO reconstruction. From 1st to 17th order (PyWENO) - + Order of the reconstruction. From 1st to 17th order (PyWENO), and from + 4 to 10 (even) in Poly ``Default = 5`` .. attribute:: time_integrator @@ -208,6 +209,9 @@ def setup(self,solution): if self.lim_type == 2: self.num_ghost = (self.reconstruction_order+1)/2 + if self.lim_type == 3: + self.num_ghost = 1 + self.reconstruction_order/2 + 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.') From e51b1d50639a09bdddd0944a8eabcc6f2558d631 Mon Sep 17 00:00:00 2001 From: damiansa Date: Tue, 20 Oct 2020 12:22:30 +0300 Subject: [PATCH 10/14] changed weno_order for reconstruction_order in function setup --- examples/advection_1d/advection_1d_nonunif.py | 2 +- examples/cubic_1d/cubic.py | 5 +++-- fvmbook/chap6/limiters.ipynb | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/advection_1d/advection_1d_nonunif.py b/examples/advection_1d/advection_1d_nonunif.py index 3df68c528..b8e014fb8 100644 --- a/examples/advection_1d/advection_1d_nonunif.py +++ b/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: diff --git a/examples/cubic_1d/cubic.py b/examples/cubic_1d/cubic.py index 53e0eee12..ac03892c1 100755 --- a/examples/cubic_1d/cubic.py +++ b/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 diff --git a/fvmbook/chap6/limiters.ipynb b/fvmbook/chap6/limiters.ipynb index 3832a9e74..472c6f9b0 100644 --- a/fvmbook/chap6/limiters.ipynb +++ b/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", From 385bc83a5e704b000f590bf2b39a886ab53bf317 Mon Sep 17 00:00:00 2001 From: damiansa Date: Tue, 20 Oct 2020 12:22:57 +0300 Subject: [PATCH 11/14] added lim_type to function setup --- examples/acoustics_1d_homogeneous/acoustics_1d.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/acoustics_1d_homogeneous/acoustics_1d.py b/examples/acoustics_1d_homogeneous/acoustics_1d.py index b06187783..c59e3e110 100644 --- a/examples/acoustics_1d_homogeneous/acoustics_1d.py +++ b/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, reconstruction_order=5, + outdir='./_output', ptwise=False, reconstruction_order=5, lim_type=2, time_integrator='SSP104', disable_output=False, output_style=1): if use_petsc: @@ -48,6 +48,7 @@ def setup(use_petsc=False, kernel_language='Fortran', solver_type='classic', solver = pyclaw.SharpClawSolver1D(riemann_solver) solver.reconstruction_order = reconstruction_order solver.time_integrator = time_integrator + solver.lim_type = lim_type if time_integrator == 'SSPLMMk3': solver.lmm_steps = 4 else: From e5ac2fe44f602f4a28af5dee36e64dd5da378275 Mon Sep 17 00:00:00 2001 From: damiansa Date: Tue, 20 Oct 2020 12:23:43 +0300 Subject: [PATCH 12/14] declared self.num_ghosts as int --- src/pyclaw/sharpclaw/solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index cfcb4e278..991ef5cac 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -207,10 +207,10 @@ 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.reconstruction_order+1)/2 + self.num_ghost = int((self.reconstruction_order+1)/2) if self.lim_type == 3: - self.num_ghost = 1 + self.reconstruction_order/2 + self.num_ghost = int(1 + self.reconstruction_order/2) 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. \ From 4de6e91af93e7aaf43653460d784fd07f2a533e4 Mon Sep 17 00:00:00 2001 From: damiansa Date: Tue, 20 Oct 2020 14:55:34 +0300 Subject: [PATCH 13/14] removed case 4 for lim_type select in reconstruct. Polynomial reconstruction can be selected using lim_type=3 --- src/pyclaw/sharpclaw/reconstruct.f90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 781b866fb..27f46905a 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -48,9 +48,6 @@ subroutine alloc_recon_workspace(maxnx,num_ghost,num_eqn,num_waves,lim_type,char case(3) allocate(uu(2,maxnx+2*num_ghost)) allocate( dq1m(maxnx+2*num_ghost)) - case(4) - allocate(uu(2,maxnx+2*num_ghost)) - allocate( dq1m(maxnx+2*num_ghost)) end select recon_alloc = .True. @@ -89,9 +86,6 @@ subroutine dealloc_recon_workspace(lim_type,char_decomp) case(3) deallocate(uu) deallocate(dq1m) - case(4) - deallocate(uu) - deallocate(dq1m) end select recon_alloc = .False. From 0f939e8c9cd655875892008be2a70cfe061360da Mon Sep 17 00:00:00 2001 From: damiansa Date: Wed, 21 Oct 2020 08:12:27 +0300 Subject: [PATCH 14/14] corrected error message in poly routine; it now describes the reconstruction_orders included' --- src/pyclaw/sharpclaw/reconstruct.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 27f46905a..c783768b5 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -849,7 +849,7 @@ subroutine poly_comp(q,ql,qr,num_eqn,maxnx,num_ghost) 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 either 4 or 6' + print *, 'ERROR: reconstruction_order for polynomial must be an odd number between 4 and 10 (inclusive)' stop end select