diff --git a/fobos b/fobos index 14c35db..9b4abc2 100644 --- a/fobos +++ b/fobos @@ -10,8 +10,7 @@ $CSHARED_INT = -cpp -c -fpic -assume realloc_lhs $LSHARED = -shared $CSTATIC_GNU = -cpp -c -frealloc-lhs $CSTATIC_INT = -cpp -c -assume realloc_lhs -$DEBUG_GNU = -Og -g3 -Warray-bounds -Wcharacter-truncation -Wline-truncation -Wimplicit-interface -Wimplicit-procedure -Wunderflow -fcheck=all -fmodule-private -ffree-line-length-132 -fimplicit-none -fbacktrace -fdump-core -finit-real=nan -; $DEBUG_GNU = -Og -g3 -Warray-bounds -Wcharacter-truncation -Wline-truncation -Wimplicit-interface -Wimplicit-procedure -Wunderflow -fcheck=all -fmodule-private -ffree-line-length-132 -fimplicit-none -fbacktrace -fdump-core -finit-real=nan -std=f2008 -fall-intrinsics +$DEBUG_GNU = -Og -g3 -Warray-bounds -Wcharacter-truncation -Wline-truncation -Wimplicit-interface -Wimplicit-procedure -Wunderflow -fcheck=all -fmodule-private -ffree-line-length-132 -fimplicit-none -fbacktrace -fdump-core -finit-real=nan -std=f2008 -fall-intrinsics $DEBUG_INT = -O0 -debug all -check all -warn all -extend-source 132 -traceback -gen-interfaces#-fpe-all=0 -fp-stack-check -fstack-protector-all -ftrapuv -no-ftz -std08 $OPTIMIZE = -O2 $EXDIRS = FLAP/src/tests/ PENF/src/tests/ pyplot-fortran/src/tests/ @@ -121,6 +120,7 @@ jobs = 2 compiler = gnu cflags = $CSHARED_GNU $DEBUG_GNU lflags = $LSHARED $DEBUG_GNU +preproc = -DDEBUG exclude_dirs = $EXDIRS mod_dir = ./mod/ obj_dir = ./obj/ @@ -134,6 +134,7 @@ jobs = 2 compiler = gnu cflags = $CSTATIC_GNU $DEBUG_GNU lflags = $DEBUG_GNU +preproc = -DDEBUG exclude_dirs = $EXDIRS mod_dir = ./mod/ obj_dir = ./obj/ @@ -173,6 +174,7 @@ jobs = 2 compiler = intel cflags = $CSHARED_INT $DEBUG_INT lflags = $LSHARED $DEBUG_INT +preproc = -DDEBUG exclude_dirs = $EXDIRS mod_dir = ./mod/ obj_dir = ./obj/ @@ -186,6 +188,7 @@ jobs = 2 compiler = intel cflags = $CSTATIC_INT $DEBUG_INT lflags = $DEBUG_INT +preproc = -DDEBUG exclude_dirs = $EXDIRS mod_dir = ./mod/ obj_dir = ./obj/ diff --git a/src/lib/abstract_objects/wenoof_alpha_object.f90 b/src/lib/abstract_objects/wenoof_alpha_object.f90 new file mode 100644 index 0000000..056d729 --- /dev/null +++ b/src/lib/abstract_objects/wenoof_alpha_object.f90 @@ -0,0 +1,40 @@ +!< Abstract alpha (non linear weights) object. +module wenoof_alpha_object +!< Abstract alpha (non linear weights) object. + +use penf, only : I_P, R_P +use wenoof_base_object +use wenoof_beta_object +use wenoof_kappa_object + +implicit none +private +public :: alpha_object +public :: alpha_object_constructor + +type, extends(base_object_constructor) :: alpha_object_constructor + !< Abstract alpha (non linear weights) object constructor. + contains +endtype alpha_object_constructor + +type, extends(base_object), abstract :: alpha_object + !< Abstract alpha (non linear weights) object. + real(R_P), allocatable :: values(:,:) !< Alpha coefficients [1:2,0:S-1]. + real(R_P), allocatable :: values_sum(:) !< Sum of alpha coefficients [1:2]. + contains + ! public deferred methods + procedure(compute_interface), pass(self), deferred :: compute !< Compute alpha. +endtype alpha_object + +abstract interface + !< Abstract interfaces of [[alpha_object]]. + pure subroutine compute_interface(self, beta, kappa) + !< Compute alpha. + import :: alpha_object, beta_object, kappa_object + class(alpha_object), intent(inout) :: self !< Alpha. + class(beta_object), intent(in) :: beta !< Beta. + class(kappa_object), intent(in) :: kappa !< Kappa. + endsubroutine compute_interface +endinterface + +endmodule wenoof_alpha_object diff --git a/src/lib/abstract_objects/wenoof_base_object.f90 b/src/lib/abstract_objects/wenoof_base_object.f90 new file mode 100644 index 0000000..00d3c5b --- /dev/null +++ b/src/lib/abstract_objects/wenoof_base_object.f90 @@ -0,0 +1,100 @@ +!< Abstract base object, the ancestor of all. +module wenoof_base_object +!< Abstract base object, the ancestor of all. +!< +!< Define a minimal, base, object that is used as ancestor of all objects, e.g. smoothness indicator, optimal weights, etc... + +use penf + +implicit none +private +public :: base_object +public :: base_object_constructor + +real(R_P), parameter :: EPS_DEF=10._R_P**(-6) !< Small epsilon to avoid division by zero, default value. + +type :: base_object_constructor + !< Abstract base object constructor. + integer(I_P) :: S=0_I_P !< Stencils dimension. + logical :: face_left=.true. !< Activate left-face interpolation computation. + logical :: face_right=.true. !< Activate right-face interpolation computation. + real(R_P) :: eps=EPS_DEF !< Small epsilon to avoid division by zero. +endtype base_object_constructor + +type, abstract :: base_object + !< Abstract base object, the ancestor of all. + !< + !< Define a minimal, base, object that is used as ancestor of all objects, e.g. smoothness indicator, optimal weights, etc... + integer(I_P) :: S=0_I_P !< Stencils dimension. + integer(I_P) :: f1=1_I_P !< Lower bound of faces index. + integer(I_P) :: f2=2_I_P !< Upper bound of faces index. + integer(I_P) :: ff=0_I_P !< Offset (step) of faces index. + real(R_P) :: eps=EPS_DEF !< Small epsilon to avoid division by zero. + contains + ! public deferred methods + procedure(create_interface), pass(self), deferred :: create !< Create object. + procedure(description_interface), pass(self), deferred :: description !< Return object string-description. + procedure(destroy_interface), pass(self), deferred :: destroy !< Destroy object. + ! public non overridable methods + procedure, pass(self), non_overridable :: create_ !< Create object. + procedure, pass(self), non_overridable :: destroy_ !< Destroy object. +endtype base_object + +abstract interface + !< Abstract interfaces of [[base_object]]. + subroutine create_interface(self, constructor) + !< Create object. + !< + !< @note Before call this method a concrete constructor must be instantiated. + import :: base_object, base_object_constructor + class(base_object), intent(inout) :: self !< Object. + class(base_object_constructor), intent(in) :: constructor !< Object constructor. + endsubroutine create_interface + + pure function description_interface(self) result(string) + !< Return object string-description. + import :: base_object + class(base_object), intent(in) :: self !< Object. + character(len=:), allocatable :: string !< String-description. + endfunction description_interface + + elemental subroutine destroy_interface(self) + !< Destroy object. + import :: base_object + class(base_object), intent(inout) :: self !< Object. + endsubroutine destroy_interface +endinterface + +contains + ! public non overridable methods + subroutine create_(self, constructor) + !< Create object. + class(base_object), intent(inout) :: self !< Object. + class(base_object_constructor), intent(in) :: constructor !< Object constructor. + + call self%destroy_ + select type(constructor) + class is(base_object_constructor) + self%S = constructor%S + if (constructor%face_left.and.constructor%face_right) then + self%f1 = 1_I_P; self%f2 = 2_I_P; self%ff = 0_I_P + elseif (constructor%face_left) then + self%f1 = 1_I_P; self%f2 = 1_I_P; self%ff = 0_I_P + elseif (constructor%face_right) then + self%f1 = 2_I_P; self%f2 = 2_I_P; self%ff = -1_I_P + endif + self%eps = constructor%eps + endselect + endsubroutine create_ + + elemental subroutine destroy_(self) + !< Destroy object. + class(base_object), intent(inout) :: self !< Object. + + self%S = 0_I_P + self%f1 = 1_I_P + self%f2 = 2_I_P + self%ff = 0_I_P + self%eps = EPS_DEF + endsubroutine destroy_ +endmodule wenoof_base_object diff --git a/src/lib/abstract_objects/wenoof_beta_object.f90 b/src/lib/abstract_objects/wenoof_beta_object.f90 new file mode 100644 index 0000000..fade5de --- /dev/null +++ b/src/lib/abstract_objects/wenoof_beta_object.f90 @@ -0,0 +1,35 @@ +!< Abstract Beta coefficients (smoothness indicators of stencil interpolations) object. +module wenoof_beta_object +!< Abstract Beta coefficients (smoothness indicators of stencil interpolations) object. + +use penf, only : I_P, R_P +use wenoof_base_object + +implicit none +private +public :: beta_object +public :: beta_object_constructor + +type, extends(base_object_constructor) :: beta_object_constructor + !< Abstract Beta coefficients object constructor. +endtype beta_object_constructor + +type, extends(base_object), abstract :: beta_object + !< Abstract Beta coefficients (smoothness indicators of stencil interpolations) object. + real(R_P), allocatable :: values(:,:) !< Beta values [1:2,0:S-1]. + contains + ! public deferred methods + procedure(compute_interface), pass(self), deferred :: compute !< Compute beta. +endtype beta_object + +abstract interface + !< Abstract interfaces of [[beta_object]]. + pure subroutine compute_interface(self, stencil) + !< Compute beta. + import :: beta_object, R_P + class(beta_object), intent(inout) :: self !< Beta. + real(R_P), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + endsubroutine compute_interface +endinterface + +endmodule wenoof_beta_object diff --git a/src/lib/abstract_objects/wenoof_interpolations_object.f90 b/src/lib/abstract_objects/wenoof_interpolations_object.f90 new file mode 100644 index 0000000..f87e706 --- /dev/null +++ b/src/lib/abstract_objects/wenoof_interpolations_object.f90 @@ -0,0 +1,35 @@ +!< Abstract interpolations object. +module wenoof_interpolations_object +!< Abstract interpolations object. + +use penf, only : I_P, R_P +use wenoof_base_object + +implicit none +private +public :: interpolations_object +public :: interpolations_object_constructor + +type, extends(base_object_constructor) :: interpolations_object_constructor + !< Abstract interpolations object constructor. +endtype interpolations_object_constructor + +type, extends(base_object), abstract :: interpolations_object + !< Abstract interpolations object. + real(R_P), allocatable :: values(:,:) !< Stencil interpolations values [1:2,0:S-1]. + contains + ! public deferred methods + procedure(compute_interface), pass(self), deferred :: compute !< Compute beta. +endtype interpolations_object + +abstract interface + !< Abstract interfaces of [[interpolations_object]]. + pure subroutine compute_interface(self, stencil) + !< Compute interpolations. + import :: interpolations_object, R_P + class(interpolations_object), intent(inout) :: self !< Interpolations. + real(R_P), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + endsubroutine compute_interface +endinterface + +endmodule wenoof_interpolations_object diff --git a/src/lib/abstract_objects/wenoof_interpolator_object.f90 b/src/lib/abstract_objects/wenoof_interpolator_object.f90 new file mode 100644 index 0000000..0e67060 --- /dev/null +++ b/src/lib/abstract_objects/wenoof_interpolator_object.f90 @@ -0,0 +1,58 @@ +!< Abstract interpolator object. +module wenoof_interpolator_object +!< Abstract interpolator object. + +use penf, only : I_P, R_P +use wenoof_base_object +use wenoof_interpolations_object +use wenoof_weights_object + +implicit none +private +public :: interpolator_object +public :: interpolator_object_constructor + +type, extends(base_object_constructor) :: interpolator_object_constructor + !< Abstract interpolator object constructor. + !< + !< @note Every concrete WENO interpolator implementations must define their own constructor type. + class(interpolations_object_constructor), allocatable :: interpolations_constructor !< Stencil interpolations constructor. + class(weights_object_constructor), allocatable :: weights_constructor !< Weights of interpolations constructor. +endtype interpolator_object_constructor + +type, extends(base_object), abstract :: interpolator_object + !< Abstract interpolator object. + !< + !< @note Do not implement any actual interpolator: provide the interface for the different interpolators implemented. + class(interpolations_object), allocatable :: interpolations !< Stencil interpolations. + class(weights_object), allocatable :: weights !< Weights of interpolations. + contains + ! public deferred methods + procedure(interpolate_debug_interface), pass(self), deferred :: interpolate_debug !< Interpolate values, debug mode. + procedure(interpolate_standard_interface), pass(self), deferred :: interpolate_standard !< Interpolate values, standard mode. + ! public methods + generic :: interpolate => interpolate_standard, interpolate_debug !< Interpolate values. +endtype interpolator_object + +abstract interface + !< Abstract interfaces of [[interpolator_object]]. + pure subroutine interpolate_debug_interface(self, stencil, interpolation, si, weights) + !< Interpolate values (providing also debug values). + import :: interpolator_object, R_P + class(interpolator_object), intent(inout) :: self !< Interpolator. + real(R_P), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. + endsubroutine interpolate_debug_interface + + pure subroutine interpolate_standard_interface(self, stencil, interpolation) + !< Interpolate values (without providing debug values). + import :: interpolator_object, R_P + class(interpolator_object), intent(inout) :: self !< Interpolator. + real(R_P), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + endsubroutine interpolate_standard_interface +endinterface + +endmodule wenoof_interpolator_object diff --git a/src/lib/abstract_objects/wenoof_kappa_object.f90 b/src/lib/abstract_objects/wenoof_kappa_object.f90 new file mode 100644 index 0000000..4f86666 --- /dev/null +++ b/src/lib/abstract_objects/wenoof_kappa_object.f90 @@ -0,0 +1,34 @@ +!< Abstract Kappa (optimal, linear weights of stencil interpolations) object. +module wenoof_kappa_object +!< Abstract Kappa (optimal, linear weights of stencil interpolations) object. + +use penf, only : I_P, R_P +use wenoof_base_object + +implicit none +private +public :: kappa_object +public :: kappa_object_constructor + +type, extends(base_object_constructor) :: kappa_object_constructor + !< Abstract kappa object constructor. +endtype kappa_object_constructor + +type, extends(base_object), abstract :: kappa_object + !< Kappa (optimal, linear weights of stencil interpolations) object. + real(R_P), allocatable :: values(:,:) !< Kappa coefficients values [1:2,0:S-1]. + contains + ! public deferred methods + procedure(compute_interface), pass(self), deferred :: compute !< Compute kappa. +endtype kappa_object + +abstract interface + !< Abstract interfaces of [[kappa_object]]. + pure subroutine compute_interface(self) + !< Compute kappa. + import :: kappa_object + class(kappa_object), intent(inout) :: self !< Kappa. + endsubroutine compute_interface +endinterface + +endmodule wenoof_kappa_object diff --git a/src/lib/abstract_objects/wenoof_weights_object.f90 b/src/lib/abstract_objects/wenoof_weights_object.f90 new file mode 100644 index 0000000..2ecae63 --- /dev/null +++ b/src/lib/abstract_objects/wenoof_weights_object.f90 @@ -0,0 +1,35 @@ +!< Abstract weights object. +module wenoof_weights_object +!< Abstract weights object. + +use penf, only : I_P, R_P +use wenoof_base_object + +implicit none +private +public :: weights_object +public :: weights_object_constructor + +type, extends(base_object_constructor) :: weights_object_constructor + !< Abstract weights object constructor. +endtype weights_object_constructor + +type, extends(base_object), abstract :: weights_object + !< Weights of stencil interpolations object. + real(R_P), allocatable :: values(:,:) !< Weights values of stencil interpolations [1:2,0:S-1]. + contains + ! deferred public methods + procedure(compute_interface), pass(self), deferred :: compute !< Compute weights. +endtype weights_object + +abstract interface + !< Abstract interfaces of [[weights_object]]. + pure subroutine compute_interface(self, stencil) + !< Compute beta. + import :: weights_object, R_P + class(weights_object), intent(inout) :: self !< Weights. + real(R_P), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + endsubroutine compute_interface +endinterface + +endmodule wenoof_weights_object diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 b/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 new file mode 100644 index 0000000..ffbfe84 --- /dev/null +++ b/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 @@ -0,0 +1,103 @@ +!< Jiang-Shu alpha (non linear weights) object. +module wenoof_alpha_rec_js +!< Jiang-Shu alpha (non linear weights) object. +!< +!< @note The provided alpha implements the alpha coefficients defined in *Efficient Implementation of Weighted ENO +!< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. + +use penf, only : I_P, R_P +use wenoof_alpha_object +use wenoof_base_object +use wenoof_beta_object +use wenoof_kappa_object + +implicit none +private +public :: alpha_rec_js +public :: alpha_rec_js_constructor +public :: create_alpha_rec_js_constructor + +type, extends(alpha_object_constructor) :: alpha_rec_js_constructor + !< Jiang-Shu alpha object constructor. +endtype alpha_rec_js_constructor + +type, extends(alpha_object) :: alpha_rec_js + !< Jiang-Shu alpha object. + !< + !< @note The provided WENO alpha implements the alpha coefficients defined in *Efficient Implementation of Weighted + !< ENO Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. + contains + ! public deferred methods + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute !< Compute alpha. + procedure, pass(self) :: description !< Return alpha string-description. + procedure, pass(self) :: destroy !< Destroy alpha. +endtype alpha_rec_js + +contains + ! public non TBP + subroutine create_alpha_rec_js_constructor(S, constructor, face_left, face_right, eps) + !< Create alpha constructor. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(alpha_object_constructor), allocatable, intent(out) :: constructor !< Constructor. + logical, intent(in), optional :: face_left !< Activate left-face interpolations. + logical, intent(in), optional :: face_right !< Activate right-face interpolations. + real(R_P), intent(in), optional :: eps !< Small epsilon to avoid division by zero. + + allocate(alpha_rec_js_constructor :: constructor) + constructor%S = S + if (present(face_left)) constructor%face_left = face_left + if (present(face_right)) constructor%face_right = face_right + if (present(eps)) constructor%eps = eps + endsubroutine create_alpha_rec_js_constructor + + ! deferred public methods + subroutine create(self, constructor) + !< Create alpha. + class(alpha_rec_js), intent(inout) :: self !< Alpha. + class(base_object_constructor), intent(in) :: constructor !< Alpha constructor. + + call self%destroy + call self%create_(constructor=constructor) + allocate(self%values(1:2, 0:self%S - 1)) + allocate(self%values_sum(1:2)) + self%values = 0._R_P + self%values_sum = 0._R_P + endsubroutine create + + pure subroutine compute(self, beta, kappa) + !< Compute alpha. + class(alpha_rec_js), intent(inout) :: self !< Alpha coefficient. + class(beta_object), intent(in) :: beta !< Beta coefficients. + class(kappa_object), intent(in) :: kappa !< Kappa coefficients. + integer(I_P) :: f, s1 !< Counters. + + self%values_sum = 0._R_P + do s1=0, self%S - 1 ! stencil loops + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + self%values(f, s1) = kappa%values(f, s1)/(self%eps + beta%values(f, s1)) ** self%S + self%values_sum(f) = self%values_sum(f) + self%values(f, s1) + enddo + enddo + endsubroutine compute + + pure function description(self) result(string) + !< Return alpha string-descripition. + class(alpha_rec_js), intent(in) :: self !< Alpha coefficient. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'alpha_rec_js%description to be implemented, do not use!' +#endif + endfunction description + + elemental subroutine destroy(self) + !< Destroy alpha. + class(alpha_rec_js), intent(inout) :: self !< Alpha. + + call self%destroy_ + if (allocated(self%values)) deallocate(self%values) + if (allocated(self%values_sum)) deallocate(self%values_sum) + endsubroutine destroy +endmodule wenoof_alpha_rec_js diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_m.F90 b/src/lib/concrete_objects/wenoof_alpha_rec_m.F90 new file mode 100644 index 0000000..3f5d5c2 --- /dev/null +++ b/src/lib/concrete_objects/wenoof_alpha_rec_m.F90 @@ -0,0 +1,110 @@ +!< Henrick alpha (non linear weights) object. +module wenoof_alpha_rec_m +!< Henrick alpha (non linear weights) object. +!< +!< @note The provided alpha implements the alpha coefficients defined in *Mapped weighted essentially non-oscillatory schemes: +!< Achieving optimal order near critical points*, Andrew K. Henrick, Tariq D. Aslam, Joseph M. Powers, JCP, +!< 2005, vol. 207, pp. 542-567, doi:10.1016/j.jcp.2005.01.023 + +use penf, only : I_P, R_P +use wenoof_alpha_object +use wenoof_alpha_rec_js +use wenoof_alpha_rec_z +use wenoof_base_object +use wenoof_beta_object +use wenoof_kappa_object + +implicit none +private +public :: alpha_rec_m +public :: alpha_rec_m_constructor + +type, extends(alpha_object_constructor) :: alpha_rec_m_constructor + !< Henrick alpha (non linear weights) object constructor. + character(len=:), allocatable :: base_type !< Base alpha coefficient type. +endtype alpha_rec_m_constructor + +type, extends(alpha_object) :: alpha_rec_m + !< Henrick alpha (non linear weights) object. + !< + !< @note The provided alpha implements the alpha coefficients defined in *Mapped weighted essentially non-oscillatory schemes: + !< Achieving optimal order near critical points*, Andrew K. Henrick, Tariq D. Aslam, Joseph M. Powers, + !< JCP, 2005, vol. 207, pp. 542-567, doi:10.1016/j.jcp.2005.01.023. + class(alpha_object), allocatable :: alpha_base !< Base alpha to be re-mapped. + contains + ! public deferred methods + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute !< Compute alpha. + procedure, pass(self) :: description !< Return alpha string-description. + procedure, pass(self) :: destroy !< Destroy alpha. +endtype alpha_rec_m + +contains + ! deferred public methods + subroutine create(self, constructor) + !< Create alpha. + class(alpha_rec_m), intent(inout) :: self !< Alpha. + class(base_object_constructor), intent(in) :: constructor !< Alpha constructor. + + call self%destroy + call self%create_(constructor=constructor) + select type(constructor) + type is(alpha_rec_m_constructor) + if (allocated(constructor%base_type)) then + select case(constructor%base_type) + case('JS') + if (allocated(self%alpha_base)) deallocate(self%alpha_base) + allocate(alpha_rec_js :: self%alpha_base) + call self%alpha_base%create(constructor=constructor) + case('Z') + if (allocated(self%alpha_base)) deallocate(self%alpha_base) + allocate(alpha_rec_z :: self%alpha_base) + call self%alpha_base%create(constructor=constructor) + endselect + endif + class default + ! @TODO add error handling + endselect + endsubroutine create + + pure subroutine compute(self, beta, kappa) + !< Compute alpha. + class(alpha_rec_m), intent(inout) :: self !< Alpha. + class(beta_object), intent(in) :: beta !< Beta. + class(kappa_object), intent(in) :: kappa !< Kappa. + integer(I_P) :: f, s1 !< Counters. + + self%values_sum = 0._R_P + call self%alpha_base%compute(beta=beta, kappa=kappa) + do s1=0, self%S - 1 ! stencil loops + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + self%values(f, s1) = & + (self%alpha_base%values(f, s1) * (kappa%values(f, s1) + kappa%values(f, s1) * kappa%values(f, s1) - & + 3._R_P * kappa%values(f, s1) * self%alpha_base%values(f, s1) + self%alpha_base%values(f, s1) * & + self%alpha_base%values(f, s1))) / & + (kappa%values(f, s1) * kappa%values(f, s1) + self%alpha_base%values(f, s1) * & + (1._R_P - 2._R_P * kappa%values(f, s1))) + self%values_sum(f) = self%values_sum(f) + self%values(f, s1) + enddo + enddo + endsubroutine compute + + pure function description(self) result(string) + !< Return alpha string-descripition. + class(alpha_rec_m), intent(in) :: self !< Alpha. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'alpha_rec_m%description to be implemented, do not use!' +#endif + endfunction description + + elemental subroutine destroy(self) + !< Destroy alpha. + class(alpha_rec_m), intent(inout) :: self !< Alpha. + + call self%destroy_ + if (allocated(self%alpha_base)) deallocate(self%alpha_base) + endsubroutine destroy +endmodule wenoof_alpha_rec_m diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_z.F90 b/src/lib/concrete_objects/wenoof_alpha_rec_z.F90 new file mode 100644 index 0000000..73c753a --- /dev/null +++ b/src/lib/concrete_objects/wenoof_alpha_rec_z.F90 @@ -0,0 +1,111 @@ +!< Borges alpha (non linear weights) object. +module wenoof_alpha_rec_z +!< Borges alpha (non linear weights) object. +!< +!< @note The provided WENO alpha implements the alpha coefficients defined in *An improved weighted essentially non-oscillatory +!< scheme for hyperbolic conservation laws*, Rafael Borges, Monique Carmona, Bruno Costa and Wai Sun Don, JCP, 2008, +!< vol. 227, pp. 3191-3211, doi: 10.1016/j.jcp.2007.11.038. + +use penf, only : I_P, R_P +use wenoof_alpha_object +use wenoof_base_object +use wenoof_beta_object +use wenoof_kappa_object + +implicit none +private +public :: alpha_rec_z +public :: alpha_rec_z_constructor + +type, extends(alpha_object_constructor) :: alpha_rec_z_constructor + !< Borges alpha (non linear weights) object constructor. +endtype alpha_rec_z_constructor + +type, extends(alpha_object) :: alpha_rec_z + !< Borges alpha (non linear weights) object. + !< + !< @note The provided alpha implements the alpha coefficients defined in *An improved weighted essentially non-oscillatory + !< scheme for hyperbolic conservation laws*, Rafael Borges, Monique Carmona, Bruno Costa and Wai Sun Don, JCP, + !< 2008, vol. 227, pp. 3191-3211, doi: 10.1016/j.jcp.2007.11.038. + contains + ! public deferred methods + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute !< Compute alpha. + procedure, pass(self) :: description !< Return alpha string-description. + procedure, pass(self) :: destroy !< Destroy alpha. +endtype alpha_rec_z +contains + ! public deferred methods + subroutine create(self, constructor) + !< Create alpha. + class(alpha_rec_z), intent(inout) :: self !< Alpha. + class(base_object_constructor), intent(in) :: constructor !< Alpha constructor. + + call self%destroy + call self%create_(constructor=constructor) + endsubroutine create + + pure subroutine compute(self, beta, kappa) + !< Compute alpha. + class(alpha_rec_z), intent(inout) :: self !< Alpha. + class(beta_object), intent(in) :: beta !< Beta. + class(kappa_object), intent(in) :: kappa !< Kappa. + integer(I_P) :: f, s1 !< Counters. + + self%values_sum = 0._R_P + do s1=0, self%S - 1 ! stencil loops + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + self%values(f, s1) = kappa%values(f, s1) * & + ((1._R_P + (tau(S=self%S, beta=beta%values) / (self%eps + beta%values(f, s1)))) ** (weno_exp(self%S))) + self%values_sum(f) = self%values_sum(f) + self%values(f, s1) + enddo + enddo + endsubroutine compute + + pure function description(self) result(string) + !< Return alpha string-descripition. + class(alpha_rec_z), intent(in) :: self !< Alpha coefficients. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'alpha_rec_z%description to be implemented, do not use!' +#endif + endfunction description + + elemental subroutine destroy(self) + !< Destroy alpha. + class(alpha_rec_z), intent(inout) :: self !< Alpha. + + call self%destroy_ + endsubroutine destroy + + ! private non TBP + pure function tau(S, beta) result(w_tau) + !< Compute the tau coefficient used in the WENO-Z alpha coefficients. + integer(I_P), intent(in) :: S !< Number of stencils used. + real(R_P), intent(in) :: beta(0:S-1) !< Smoothness indicators. + real(R_P) :: w_tau !< Tau coefficient. + + w_tau = abs(beta(0) - & + (1_I_P - weno_odd(S)) * beta(1) - & + (1_I_P - weno_odd(S)) * beta(S-2) + & + (1_I_P - 2_I_P * weno_odd(S)) * beta(S-1)) + endfunction tau + + pure function weno_exp(S) result(w_exp) + !< Compute the exponent used in the alpha function. + integer(I_P), intent(in) :: S !< Number of stencils used. + integer(I_P) :: w_exp !< Exponent used in the alpha function. + + w_exp = int(S, I_P) + endfunction weno_exp + + pure function weno_odd(S) result(w_odd) + !< Compute the distinguisher between odd and even number of stencils. + integer(I_P), intent(in) :: S !< Number of stencils used. + integer(I_P) :: w_odd !< Distinguishing between odd and even number of stencils. + + w_odd = int(mod(S, 2_I_P), I_P) + endfunction weno_odd +endmodule wenoof_alpha_rec_z diff --git a/src/lib/wenoof_smoothness_indicators_js.f90 b/src/lib/concrete_objects/wenoof_beta_rec_js.F90 similarity index 97% rename from src/lib/wenoof_smoothness_indicators_js.f90 rename to src/lib/concrete_objects/wenoof_beta_rec_js.F90 index f85cda4..1ad9b6f 100644 --- a/src/lib/wenoof_smoothness_indicators_js.f90 +++ b/src/lib/concrete_objects/wenoof_beta_rec_js.F90 @@ -1,113 +1,71 @@ -!< Jiang-Shu and Gerolymos-Senechal-Vallet smoothness indicators object. -module wenoof_smoothness_indicators_js -!< Jiang-Shu and Gerolymos-Senechal-Vallet smoothness indicators object. +!< Jiang-Shu and Gerolymos-Senechal-Vallet Beta coefficients (smoothness indicators of stencil interpolations) object. +module wenoof_beta_rec_js +!< Jiang-Shu and Gerolymos-Senechal-Vallet Beta coefficients (smoothness indicators of stencil interpolations) object. !< -!< @note The provided WENO optimal weights implements the smoothness indicators defined in *Efficient Implementation +!< @note The provided beta object implements the smoothness indicators defined in *Efficient Implementation !< of Weighted ENO Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, !< doi:10.1016/j.jcp.2009.07.039 use penf, only : I_P, R_P use wenoof_base_object -use wenoof_smoothness_indicators +use wenoof_beta_object implicit none private -public :: smoothness_indicators_js -public :: smoothness_indicators_js_constructor -public :: create_smoothness_indicators_js_constructor +public :: beta_rec_js +public :: beta_rec_js_constructor +public :: create_beta_rec_js_constructor -type, extends(smoothness_indicators_constructor) :: smoothness_indicators_js_constructor - !< Jiang-Shu and Gerolymos-Senechal-Vallet smoothness indicators object constructor. -endtype smoothness_indicators_js_constructor +type, extends(beta_object_constructor) :: beta_rec_js_constructor + !< Jiang-Shu and Gerolymos-Senechal-Vallet beta object constructor. +endtype beta_rec_js_constructor -type, extends(smoothness_indicators) :: smoothness_indicators_js - !< Jiang-Shu and Gerolymos-Senechal-Vallet smoothness indicators object. +type, extends(beta_object) :: beta_rec_js + !< Jiang-Shu and Gerolymos-Senechal-Vallet Beta coefficients (smoothness indicators of stencil interpolations) object. !< - !< @note The provided WENO optimal weights implements the optimal weights defined in *Efficient Implementation of Weighted ENO + !< @note The provided beta object implements the smoothness indicators defined in *Efficient Implementation of Weighted ENO !< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and - !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Sénéchal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, + !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, !< doi:10.1016/j.jcp.2009.07.039 private - real(R_P), allocatable :: coef(:,:,:) !< Smoothness indicators coefficients [1:2,0:S-1,0:S-1]. + real(R_P), allocatable :: coef(:,:,:) !< Beta coefficients [1:2,0:S-1,0:S-1]. contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute smoothness indicators. - procedure, nopass :: description !< Return smoothness indicators string-description. - ! overridden public methods - procedure, pass(self) :: create !< Create smoothness indicators. - procedure, pass(self) :: destroy !< Destroy smoothness indicators. -endtype smoothness_indicators_js + ! public deferred methods + procedure, pass(self) :: create !< Create beta. + procedure, pass(self) :: compute !< Compute beta. + procedure, pass(self) :: description !< Return beta string-description. + procedure, pass(self) :: destroy !< Destroy beta. +endtype beta_rec_js contains - ! public non TBP - subroutine create_smoothness_indicators_js_constructor(S, constructor) - !< Create smoothness indicators constructor. - integer(I_P), intent(in) :: S !< Stencils dimension. - class(smoothness_indicators_constructor), allocatable, intent(out) :: constructor !< Smoothness indicators constructor. - - allocate(smoothness_indicators_js_constructor :: constructor) + ! public non TBP procedures + subroutine create_beta_rec_js_constructor(S, constructor, face_left, face_right) + !< Create beta constructor. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(beta_object_constructor), allocatable, intent(out) :: constructor !< Constructor. + logical, intent(in), optional :: face_left !< Activate left-face interpolations. + logical, intent(in), optional :: face_right !< Activate right-face interpolations. + + allocate(beta_rec_js_constructor :: constructor) constructor%S = S - endsubroutine create_smoothness_indicators_js_constructor - - ! deferred public methods - pure subroutine compute(self, S, stencil, f1, f2, ff) - !< Compute smoothness indicators. - class(smoothness_indicators_js), intent(inout) :: self !< Smoothness indicator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed. - integer(I_P) :: s1, s2, s3, f !< Counters - - do s1=0, S - 1 ! stencils loop - do f=f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - self%si(f, s1) = 0._R_P - do s2=0, S - 1 - do s3=0, S - 1 - self%si(f, s1) = self%si(f, s1) + self%coef(s3, s2, s1) * stencil(f + ff, s1 - s3) * stencil(f + ff, s1 - s2) - enddo - enddo - enddo - enddo - endsubroutine compute - - pure function description() result(string) - !< Return smoothness indicators string-description. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - - string = 'WENO smoothness indicators'//nl - string = string//' Based on the work by Jiang and Shu "Efficient Implementation of Weighted ENO Schemes", see '// & - 'JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and'//nl - string = string//' on the work by Gerolimos, Sénéchal and Vallet "Very-High-Order WENO Schemes", see '// & - 'JCP, 2009, vol. 228, pp. 8481-8524, doi:10.1016/j.jcp.2009.07.039'//nl - string = string//' The "compute" method has the following public API'//nl - string = string//' compute(S, stencil, f1, f2, ff)'//nl - string = string//' where:'//nl - string = string//' S: integer(I_P), intent(in), the number of the stencils used'//nl - string = string//' stencil: real(R_P), intent(IN), the stencil used for the interpolation [1:2, 1-S:-1+S]'//nl - string = string//' f1, f2: integer(I_P), intent(in), the faces to be computed (1 => left interface, 2 => right interface)'//nl - string = string//' ff: integer(I_P), intent(in), the parameter for the stencil value choice' - endfunction description + if (present(face_left)) constructor%face_left = face_left + if (present(face_right)) constructor%face_right = face_right + endsubroutine create_beta_rec_js_constructor - ! overridden public methods - pure subroutine create(self, constructor) - !< Create smoothness indicators. - class(smoothness_indicators_js), intent(inout) :: self !< Smoothness indicators. - class(base_object_constructor), intent(in) :: constructor !< Smoothness indicators constructor. - integer(I_P) :: S !< Stencils dimension. + ! public deferred methods + subroutine create(self, constructor) + !< Create beta. + class(beta_rec_js), intent(inout) :: self !< Beta. + class(base_object_constructor), intent(in) :: constructor !< Beta constructor. call self%destroy - call self%smoothness_indicators%create(constructor=constructor) - select type(constructor) - class is(smoothness_indicators_js_constructor) - S = constructor%S - allocate(self%coef(0:S - 1, 0:S - 1, 0:S - 1)) - class default - ! @TODO add error handling - endselect + call self%create_(constructor=constructor) + allocate(self%values(1:2, 0:self%S - 1)) + self%values = 0._R_P + allocate(self%coef(0:self%S - 1, 0:self%S - 1, 0:self%S - 1)) associate(c => self%coef) - select case(S) + select case(self%S) case(2) ! 3rd order ! stencil 0 ! i*i ; (i-1)*i @@ -2425,11 +2383,42 @@ pure subroutine create(self, constructor) endassociate endsubroutine create + pure subroutine compute(self, stencil) + !< Compute beta. + class(beta_rec_js), intent(inout) :: self !< Beta. + real(R_P), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + integer(I_P) :: s1, s2, s3, f !< Counters. + + do s1=0, self%S - 1 ! stencils loop + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + self%values(f, s1) = 0._R_P + do s2=0, self%S - 1 + do s3=0, self%S - 1 + self%values(f, s1) = self%values(f, s1) + & + self%coef(s3, s2, s1) * stencil(f + self%ff, s1 - s3) * stencil(f + self%ff, s1 - s2) + enddo + enddo + enddo + enddo + endsubroutine compute + + pure function description(self) result(string) + !< Return beta string-description. + class(beta_rec_js), intent(in) :: self !< Beta. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'beta_rec_js%description to be implemented, do not use!' +#endif + endfunction description + elemental subroutine destroy(self) - !< Destroy smoothness indicators. - class(smoothness_indicators_js), intent(inout) :: self !< Smoothenss indicators. + !< Destroy beta. + class(beta_rec_js), intent(inout) :: self !< Beta. - call self%smoothness_indicators%destroy + call self%destroy_ + if (allocated(self%values)) deallocate(self%values) if (allocated(self%coef)) deallocate(self%coef) endsubroutine destroy -endmodule wenoof_smoothness_indicators_js +endmodule wenoof_beta_rec_js diff --git a/src/lib/wenoof_polynomials_js.f90 b/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 similarity index 85% rename from src/lib/wenoof_polynomials_js.f90 rename to src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 index d14e54d..1ebbc51 100644 --- a/src/lib/wenoof_polynomials_js.f90 +++ b/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 @@ -1,111 +1,71 @@ -!< Jiang-Shu (Lagrange) polynomials object. -module wenoof_polynomials_js -!< Jiang-Shu (Lagrange) polynomials object. +!< Jiang-Shu (Lagrange) interpolations object for derivative reconstruction. +module wenoof_interpolations_rec_js +!< Jiang-Shu (Lagrange) interpolations object for derivative reconstruction. !< -!< @note The provided polynomials implement the Lagrange polynomials defined in *Efficient Implementation +!< @note The provided interpolations implement the Lagrange interpolations defined in *Efficient Implementation !< of Weighted ENO Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, !< doi:10.1016/j.jcp.2009.07.039 use penf, only : I_P, R_P use wenoof_base_object -use wenoof_polynomials +use wenoof_interpolations_object implicit none private -public :: polynomials_js -public :: polynomials_js_constructor -public :: create_polynomials_js_constructor +public :: interpolations_rec_js +public :: interpolations_rec_js_constructor +public :: create_interpolations_rec_js_constructor -type, extends(polynomials_constructor) :: polynomials_js_constructor - !< Jiang-Shu (Lagrange) polynomials object constructor. -endtype polynomials_js_constructor +type, extends(interpolations_object_constructor) :: interpolations_rec_js_constructor + !< Jiang-Shu (Lagrange) interpolations object for derivative reconstruction constructor. +endtype interpolations_rec_js_constructor -type, extends(polynomials) :: polynomials_js - !< Jiang-Shu (Lagrange) polynomials object. +type, extends(interpolations_object) :: interpolations_rec_js + !< Jiang-Shu (Lagrange) interpolations object for derivative reconstruction. !< - !< @note The provided polynomials implement the Lagrange polynomials defined in *Efficient Implementation + !< @note The provided interpolations implement the Lagrange interpolations defined in *Efficient Implementation !< of Weighted ENO Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, !< doi:10.1016/j.jcp.2009.07.039 private real(R_P), allocatable :: coef(:,:,:) !< Polynomial coefficients [1:2,0:S-1,0:S-1]. contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute weights. - procedure, nopass :: description !< Return weights string-description. - ! overridden public methods - procedure, pass(self) :: create !< Create weights. - procedure, pass(self) :: destroy !< Destroy weights. -endtype polynomials_js + ! public deferred methods + procedure, pass(self) :: create !< Create interpolations. + procedure, pass(self) :: compute !< Compute interpolations. + procedure, pass(self) :: description !< Return interpolations string-description. + procedure, pass(self) :: destroy !< Destroy interpolations. +endtype interpolations_rec_js contains - ! public non TBP - subroutine create_polynomials_js_constructor(S, constructor) - !< Create polynomials constructor. - integer(I_P), intent(in) :: S !< Stencils dimension. - class(polynomials_constructor), allocatable, intent(out) :: constructor !< Polynomials constructor. + ! public non TBP procedures + subroutine create_interpolations_rec_js_constructor(S, constructor, face_left, face_right) + !< Create interpolations constructor. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(interpolations_object_constructor), allocatable, intent(out) :: constructor !< Constructor. + logical, intent(in), optional :: face_left !< Activate left-face interpolations. + logical, intent(in), optional :: face_right !< Activate right-face interpolations. - allocate(polynomials_js_constructor :: constructor) + allocate(interpolations_rec_js_constructor :: constructor) constructor%S = S - endsubroutine create_polynomials_js_constructor + if (present(face_left)) constructor%face_left = face_left + if (present(face_right)) constructor%face_right = face_right + endsubroutine create_interpolations_rec_js_constructor - ! deferred public methods - pure subroutine compute(self, S, stencil, f1, f2, ff) - !< Compute polynomials. - class(polynomials_js), intent(inout) :: self !< WENO polynomial. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed. - integer(I_P) :: s1, s2, f !< Counters - - self%poly = 0._R_P - do s1 = 0, S - 1 ! stencils loop - do s2 = 0, S - 1 ! values loop - do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - self%poly(f, s1) = self%poly(f, s1) + self%coef(f, s2, s1) * stencil(f + ff, -s2 + s1) - enddo - enddo - enddo - endsubroutine compute - - pure function description() result(string) - !< Return polynomials string-description. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - - string = 'WENO polynomial'//nl - string = string//' Based on the work by Jiang and Shu "Efficient Implementation of Weighted ENO Schemes", see '// & - 'JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and'//nl - string = string//' on the work by Gerolimos, Sénéchal and Vallet "Very-High-Order WENO Schemes", see '// & - 'JCP, 2009, vol. 228, pp. 8481-8524, doi:10.1016/j.jcp.2009.07.039'//nl - string = string//' The "compute" method has the following public API'//nl - string = string//' compute(S, stencil, f1, f2, ff)'//nl - string = string//' where:'//nl - string = string//' S: integer(I_P), intent(in), the number of the stencils used'//nl - string = string//' stencil: real(R_P), intent(IN), the stencil used for the interpolation [1:2, 1-S:-1+S]'//nl - string = string//' f1, f2: integer(I_P), intent(in), the faces to be computed (1 => left interface, 2 => right interface)'//nl - string = string//' ff: integer(I_P), intent(in), the parameter for the stencil value choice' - endfunction description - - ! overridden public methods - pure subroutine create(self, constructor) - !< Create coefficients. - class(polynomials_js), intent(inout) :: self !< Polynomials. - class(base_object_constructor), intent(in) :: constructor !< Polynomials constructor. - integer(I_P) :: S !< Stencils dimension. + ! public deferred methods + subroutine create(self, constructor) + !< Create interpolations. + class(interpolations_rec_js), intent(inout) :: self !< Interpolations. + class(base_object_constructor), intent(in) :: constructor !< Interpolations constructor. call self%destroy - call self%polynomials%create(constructor=constructor) - select type(constructor) - class is(polynomials_js_constructor) - S = constructor%S - allocate(self%coef(1:2, 0:S - 1, 0:S - 1)) - class default - ! @TODO add error handling - endselect + call self%create_(constructor=constructor) + allocate(self%values(1:2, 0:self%S - 1)) + self%values = 0._R_P + allocate(self%coef(1:2, 0:self%S - 1, 0:self%S - 1)) associate(c => self%coef) - select case(S) + select case(self%S) case(2) ! 3rd order ! 1 => left interface (i-1/2) ! cell 0 ; cell 1 @@ -372,11 +332,41 @@ pure subroutine create(self, constructor) endassociate endsubroutine create + pure subroutine compute(self, stencil) + !< Compute interpolations. + class(interpolations_rec_js), intent(inout) :: self !< Interpolations. + real(R_P), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + integer(I_P) :: s1 !< Counter. + integer(I_P) :: s2 !< Counter. + integer(I_P) :: f !< Counter. + + self%values = 0._R_P + do s1=0, self%S - 1 ! stencils loop + do s2=0, self%S - 1 ! values loop + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + self%values(f, s1) = self%values(f, s1) + self%coef(f, s2, s1) * stencil(f + self%ff, -s2 + s1) + enddo + enddo + enddo + endsubroutine compute + + pure function description(self) result(string) + !< Return interpolations string-description. + class(interpolations_rec_js), intent(in) :: self !< Interpolations. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'interpolations_rec_js%description to be implemented, do not use!' +#endif + endfunction description + elemental subroutine destroy(self) - !< Destroy polynomials. - class(polynomials_js), intent(inout) :: self !< WENO polynomials. + !< Destroy interpolations. + class(interpolations_rec_js), intent(inout) :: self !< Interpolations. - call self%polynomials%destroy + call self%destroy_ + if (allocated(self%values)) deallocate(self%values) if (allocated(self%coef)) deallocate(self%coef) endsubroutine destroy -endmodule wenoof_polynomials_js +endmodule wenoof_interpolations_rec_js diff --git a/src/lib/concrete_objects/wenoof_interpolator_js.F90 b/src/lib/concrete_objects/wenoof_interpolator_js.F90 new file mode 100644 index 0000000..504989c --- /dev/null +++ b/src/lib/concrete_objects/wenoof_interpolator_js.F90 @@ -0,0 +1,67 @@ +!< Jiang-Shu (upwind) interpolator object. +module wenoof_interpolator_js +!< Jiang-Shu (upwind) interpolator object. + +use, intrinsic :: iso_fortran_env, only : stderr=>error_unit +use penf, only : I_P, R_P, str +use wenoof_interpolator_object + +implicit none +private +public :: interpolator_js +public :: interpolator_js_constructor + +type, extends(interpolator_object_constructor) :: interpolator_js_constructor + !< Jiang-Shu (upwind) interpolator object constructor. +endtype interpolator_js_constructor + +type, extends(interpolator_object) :: interpolator_js + !< Jiang-Shu (upwind) interpolator object. + contains + ! public deferred methods + procedure, pass(self) :: description !< Return interpolator string-description. + procedure, pass(self) :: interpolate_standard !< Interpolate values (without providing debug values). + procedure, pass(self) :: interpolate_debug !< Interpolate values (providing also debug values). +endtype interpolator_js + +contains + ! public deferred methods + pure function description(self) result(string) + !< Return interpolator string-descripition. + class(interpolator_js), intent(in) :: self !< Interpolator. + character(len=:), allocatable :: string !< String-description. + character(len=1), parameter :: nl=new_line('a') !< New line character. + character(len=:), allocatable :: dummy_string !< Dummy string. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'interpolator_js to be implemented, do not use!' +#endif + endfunction description + + pure subroutine interpolate_standard(self, stencil, interpolation) + !< Interpolate values (without providing debug values). + class(interpolator_js), intent(inout) :: self !< Interpolator. + real(R_P), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'interpolator_js to be implemented, do not use!' +#endif + endsubroutine interpolate_standard + + pure subroutine interpolate_debug(self, stencil, interpolation, si, weights) + !< Interpolate values (providing also debug values). + class(interpolator_js), intent(inout) :: self !< Interpolator. + real(R_P), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'interpolator_js to be implemented, do not use!' +#endif + endsubroutine interpolate_debug +endmodule wenoof_interpolator_js diff --git a/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 b/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 new file mode 100644 index 0000000..1cbedab --- /dev/null +++ b/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 @@ -0,0 +1,205 @@ +!< Jiang-Shu and Gerolymos-Senechal-Vallet kappa coefficients for reconstruction. +module wenoof_kappa_rec_js +!< Jiang-Shu and Gerolymos-Senechal-Vallet kappa coefficients for reconstruction. +!< +!< @note The provided WENO kappa implements the linear weights defined in *Efficient Implementation of Weighted ENO +!< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and +!< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, +!< doi:10.1016/j.jcp.2009.07.039 + +use penf, only : I_P, R_P +use wenoof_base_object +use wenoof_kappa_object + +implicit none +private +public :: kappa_rec_js +public :: kappa_rec_js_constructor +public :: create_kappa_rec_js_constructor + +type, extends(kappa_object_constructor) :: kappa_rec_js_constructor + !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal kappa object constructor. +endtype kappa_rec_js_constructor + +type, extends(kappa_object):: kappa_rec_js + !< Jiang-Shu and Gerolymos-Senechal-Vallet kappa object. + !< + !< @note The provided WENO kappa implements the weights defined in *Efficient Implementation of Weighted ENO + !< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and + !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, + !< doi:10.1016/j.jcp.2009.07.039 + contains + ! public deferred methods + procedure, pass(self) :: create !< Create kappa. + procedure, pass(self) :: compute !< Compute kappa. + procedure, pass(self) :: description !< Return kappa string-description. + procedure, pass(self) :: destroy !< Destroy kappa. +endtype kappa_rec_js + +contains + ! public non TBP procedures + subroutine create_kappa_rec_js_constructor(S, constructor) + !< Create kappa constructor. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(kappa_object_constructor), allocatable, intent(out) :: constructor !< Constructor. + + allocate(kappa_rec_js_constructor :: constructor) + constructor%S = S + endsubroutine create_kappa_rec_js_constructor + + ! deferred public methods + subroutine create(self, constructor) + !< Create kappa. + !< + !< @note The kappa coefficients are also computed, they being constants. + class(kappa_rec_js), intent(inout) :: self !< Kappa. + class(base_object_constructor), intent(in) :: constructor !< Kappa constructor. + + call self%destroy + call self%create_(constructor=constructor) + allocate(self%values(1:2, 0:self%S - 1)) + self%values = 0._R_P + call self%compute + endsubroutine create + + pure subroutine compute(self) + !< Compute kappa. + class(kappa_rec_js), intent(inout) :: self !< Kappa. + + associate(val => self%values) + select case(self%S) + case(2) ! 3rd order + ! 1 => left interface (i-1/2) + val(1, 0) = 2._R_P/3._R_P ! stencil 0 + val(1, 1) = 1._R_P/3._R_P ! stencil 1 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/3._R_P ! stencil 0 + val(2, 1) = 2._R_P/3._R_P ! stencil 1 + case(3) ! 5th order + ! 1 => left interface (i-1/2) + val(1, 0) = 0.3_R_P ! stencil 0 + val(1, 1) = 0.6_R_P ! stencil 1 + val(1, 2) = 0.1_R_P ! stencil 2 + ! 2 => right interface (i+1/2) + val(2, 0) = 0.1_R_P ! stencil 0 + val(2, 1) = 0.6_R_P ! stencil 1 + val(2, 2) = 0.3_R_P ! stencil 2 + case(4) ! 7th order + ! 1 => left interface (i-1/2) + val(1, 0) = 4._R_P/35._R_P ! stencil 0 + val(1, 1) = 18._R_P/35._R_P ! stencil 1 + val(1, 2) = 12._R_P/35._R_P ! stencil 2 + val(1, 3) = 1._R_P/35._R_P ! stencil 3 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/35._R_P ! stencil 0 + val(2, 1) = 12._R_P/35._R_P ! stencil 1 + val(2, 2) = 18._R_P/35._R_P ! stencil 2 + val(2, 3) = 4._R_P/35._R_P ! stencil 3 + case(5) ! 9th order + ! 1 => left interface (i-1/2) + val(1, 0) = 5._R_P/126._R_P ! stencil 0 + val(1, 1) = 20._R_P/63._R_P ! stencil 1 + val(1, 2) = 10._R_P/21._R_P ! stencil 2 + val(1, 3) = 10._R_P/63._R_P ! stencil 3 + val(1, 4) = 1._R_P/126._R_P ! stencil 4 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/126._R_P ! stencil 0 + val(2, 1) = 10._R_P/63._R_P ! stencil 1 + val(2, 2) = 10._R_P/21._R_P ! stencil 2 + val(2, 3) = 20._R_P/63._R_P ! stencil 3 + val(2, 4) = 5._R_P/126._R_P ! stencil 4 + case(6) ! 11th order + ! 1 => left interface (i-1/2) + val(1, 0) = 1._R_P/77._R_P ! stencil 0 + val(1, 1) = 25._R_P/154._R_P ! stencil 1 + val(1, 2) = 100._R_P/231._R_P ! stencil 2 + val(1, 3) = 25._R_P/77._R_P ! stencil 3 + val(1, 4) = 5._R_P/77._R_P ! stencil 4 + val(1, 5) = 1._R_P/462._R_P ! stencil 5 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/462._R_P ! stencil 0 + val(2, 1) = 5._R_P/77._R_P ! stencil 1 + val(2, 2) = 25._R_P/77._R_P ! stencil 2 + val(2, 3) = 100._R_P/231._R_P ! stencil 3 + val(2, 4) = 25._R_P/154._R_P ! stencil 4 + val(2, 5) = 1._R_P/77._R_P ! stencil 5 + case(7) ! 13th order + ! 1 => left interface (i-1/2) + val(1, 0) = 7._R_P/1716._R_P ! stencil 0 + val(1, 1) = 21._R_P/286._R_P ! stencil 1 + val(1, 2) = 175._R_P/572._R_P ! stencil 2 + val(1, 3) = 175._R_P/429._R_P ! stencil 3 + val(1, 4) = 105._R_P/572._R_P ! stencil 4 + val(1, 5) = 7._R_P/286._R_P ! stencil 5 + val(1, 6) = 1._R_P/1716._R_P ! stencil 6 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/1716._R_P ! stencil 0 + val(2, 1) = 7._R_P/286._R_P ! stencil 1 + val(2, 2) = 105._R_P/572._R_P ! stencil 2 + val(2, 3) = 175._R_P/429._R_P ! stencil 3 + val(2, 4) = 175._R_P/572._R_P ! stencil 4 + val(2, 5) = 21._R_P/286._R_P ! stencil 5 + val(2, 6) = 7._R_P/1716._R_P ! stencil 6 + case(8) ! 15th order + ! 1 => left interface (i-1/2) + val(1, 0) = 8._R_P/6435._R_P ! stencil 0 + val(1, 1) = 196._R_P/6435._R_P ! stencil 1 + val(1, 2) = 392._R_P/2145._R_P ! stencil 2 + val(1, 3) = 490._R_P/1287._R_P ! stencil 3 + val(1, 4) = 392._R_P/1287._R_P ! stencil 4 + val(1, 5) = 196._R_P/2145._R_P ! stencil 5 + val(1, 6) = 56._R_P/6435._R_P ! stencil 6 + val(1, 7) = 1._R_P/6435._R_P ! stencil 7 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/6435._R_P ! stencil 0 + val(2, 1) = 56._R_P/6435._R_P ! stencil 1 + val(2, 2) = 196._R_P/2145._R_P ! stencil 2 + val(2, 3) = 392._R_P/1287._R_P ! stencil 3 + val(2, 4) = 490._R_P/1287._R_P ! stencil 4 + val(2, 5) = 392._R_P/2145._R_P ! stencil 5 + val(2, 6) = 196._R_P/6435._R_P ! stencil 6 + val(2, 7) = 8._R_P/6435._R_P ! stencil 7 + case(9) ! 17th order + ! 1 => left interface (i-1/2) + val(1, 0) = 9._R_P/24310._R_P ! stencil 0 + val(1, 1) = 144._R_P/12155._R_P ! stencil 1 + val(1, 2) = 1176._R_P/12155._R_P ! stencil 2 + val(1, 3) = 3528._R_P/12155._R_P ! stencil 3 + val(1, 4) = 882._R_P/2431._R_P ! stencil 4 + val(1, 5) = 2352._R_P/12155._R_P ! stencil 5 + val(1, 6) = 504._R_P/12155._R_P ! stencil 6 + val(1, 7) = 36._R_P/12155._R_P ! stencil 7 + val(1, 8) = 1._R_P/24310._R_P ! stencil 8 + ! 2 => right interface (i+1/2) + val(2, 0) = 1._R_P/24310._R_P ! stencil 0 + val(2, 1) = 36._R_P/12155._R_P ! stencil 1 + val(2, 2) = 504._R_P/12155._R_P ! stencil 2 + val(2, 3) = 2352._R_P/12155._R_P ! stencil 3 + val(2, 4) = 882._R_P/2431._R_P ! stencil 4 + val(2, 5) = 3528._R_P/12155._R_P ! stencil 5 + val(2, 6) = 1176._R_P/12155._R_P ! stencil 6 + val(2, 7) = 144._R_P/12155._R_P ! stencil 7 + val(2, 8) = 9._R_P/24310._R_P ! stencil 8 + endselect + endassociate + endsubroutine compute + + pure function description(self) result(string) + !< Return string-description of kappa. + class(kappa_rec_js), intent(in) :: self !< Kappa. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'kappa_rec_js%description to be implemented, do not use!' +#endif + endfunction description + + elemental subroutine destroy(self) + !< Destroy kappa. + class(kappa_rec_js), intent(inout) :: self !< Kappa. + + call self%destroy_ + if (allocated(self%values)) deallocate(self%values) + endsubroutine destroy +endmodule wenoof_kappa_rec_js diff --git a/src/lib/concrete_objects/wenoof_reconstructor_js.F90 b/src/lib/concrete_objects/wenoof_reconstructor_js.F90 new file mode 100644 index 0000000..3eeb58e --- /dev/null +++ b/src/lib/concrete_objects/wenoof_reconstructor_js.F90 @@ -0,0 +1,133 @@ +!< Jiang-Shu (upwind) reconstructor object. +module wenoof_reconstructor_js +!< Jiang-Shu (upwind) reconstructor object. + +use, intrinsic :: iso_fortran_env, only : stderr=>error_unit +use penf, only : I_P, R_P, str +use wenoof_base_object +use wenoof_interpolations_object +use wenoof_interpolations_rec_js +use wenoof_interpolator_object +use wenoof_objects_factory +use wenoof_weights_object +use wenoof_weights_js + +implicit none +private +public :: reconstructor_js +public :: reconstructor_js_constructor +public :: create_reconstructor_js_constructor + +type, extends(interpolator_object_constructor) :: reconstructor_js_constructor + !< Jiang-Shu (upwind) reconstructor object constructor. +endtype reconstructor_js_constructor + +type, extends(interpolator_object) :: reconstructor_js + !< Jiang-Shu (upwind) reconstructor object. + !< + !< @note Provide the *Efficient Implementation of Weighted ENO Schemes*, + !< Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. + !< + !< @note The supported accuracy formal order are: 3rd, 5th, 7th, 9th, 11th, 13th, 15th, 17th corresponding to use 2, 3, 4, 5, 6, + !< 7, 8, 9 stencils composed of 2, 3, 4, 5, 6, 7, 8, 9 values, respectively. + contains + ! public deferred methods + procedure, pass(self) :: create !< Create reconstructor. + procedure, pass(self) :: description !< Return reconstructor string-description. + procedure, pass(self) :: destroy !< Destroy reconstructor. + procedure, pass(self) :: interpolate_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: interpolate_standard !< Interpolate values (without providing debug values). +endtype reconstructor_js + +contains + ! public non TBP + subroutine create_reconstructor_js_constructor(S, interpolations_constructor, weights_constructor, constructor, & + face_left, face_right, eps) + !< Create weights constructor. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(interpolations_object_constructor), intent(in) :: interpolations_constructor !< Interpolations constructor. + class(weights_object_constructor), intent(in) :: weights_constructor !< Weights constructor. + class(interpolator_object_constructor), intent(out), allocatable :: constructor !< Constructor. + logical, intent(in), optional :: face_left !< Activate left-face interp. + logical, intent(in), optional :: face_right !< Activate right-face interp. + real(R_P), intent(in), optional :: eps !< Small eps to avoid zero-div. + + allocate(reconstructor_js_constructor :: constructor) + constructor%S = S + if (present(face_left)) constructor%face_left = face_left + if (present(face_right)) constructor%face_right = face_right + if (present(eps)) constructor%eps = eps + select type(constructor) + type is(reconstructor_js_constructor) + allocate(constructor%interpolations_constructor, source=interpolations_constructor) + allocate(constructor%weights_constructor, source=weights_constructor) + endselect + endsubroutine create_reconstructor_js_constructor + + ! public deferred methods + subroutine create(self, constructor) + !< Create interpolator. + class(reconstructor_js), intent(inout) :: self !< Interpolator. + class(base_object_constructor), intent(in) :: constructor !< Constructor. + type(objects_factory) :: factory !< Objects factory. + + call self%destroy + call self%create_(constructor=constructor) + select type(constructor) + type is(reconstructor_js_constructor) + call factory%create(constructor=constructor%interpolations_constructor, object=self%interpolations) + call factory%create(constructor=constructor%weights_constructor, object=self%weights) + endselect + endsubroutine create + + pure function description(self) result(string) + !< Return reconstructor string-descripition. + class(reconstructor_js), intent(in) :: self !< Reconstructor. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'reconstructor_js%description to be implemented, do not use!' +#endif + endfunction description + + elemental subroutine destroy(self) + !< Destroy reconstructor. + class(reconstructor_js), intent(inout) :: self !< Reconstructor. + + call self%destroy_ + if (allocated(self%interpolations)) deallocate(self%interpolations) + if (allocated(self%weights)) deallocate(self%weights) + endsubroutine destroy + + pure subroutine interpolate_debug(self, stencil, interpolation, si, weights) + !< Interpolate values (providing also debug values). + class(reconstructor_js), intent(inout) :: self !< Reconstructor. + real(R_P), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. + + call self%interpolate_standard(stencil=stencil, interpolation=interpolation) + si = 0._R_P ! @TODO implement beta extraction + weights = self%weights%values + endsubroutine interpolate_debug + + pure subroutine interpolate_standard(self, stencil, interpolation) + !< Interpolate values (without providing debug values). + class(reconstructor_js), intent(inout) :: self !< Reconstructor. + real(R_P), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + integer(I_P) :: f, s !< Counters. + + call self%interpolations%compute(stencil=stencil) + call self%weights%compute(stencil=stencil) + interpolation = 0._R_P + do s=0, self%S - 1 ! stencils loop + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + interpolation(f + self%ff) = interpolation(f + self%ff) + & + self%weights%values(f + self%ff, s) * self%interpolations%values(f, s) + enddo + enddo + endsubroutine interpolate_standard +endmodule wenoof_reconstructor_js diff --git a/src/lib/concrete_objects/wenoof_weights_js.F90 b/src/lib/concrete_objects/wenoof_weights_js.F90 new file mode 100644 index 0000000..95fc33f --- /dev/null +++ b/src/lib/concrete_objects/wenoof_weights_js.F90 @@ -0,0 +1,159 @@ +!< Jiang-Shu and Gerolymos-Senechal-Vallet weights. +module wenoof_weights_js +!< Jiang-Shu and Gerolymos-Senechal-Vallet weights. +!< +!< @note The provided WENO weights implements the weights defined in *Efficient Implementation of Weighted ENO +!< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and +!< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, +!< doi:10.1016/j.jcp.2009.07.039 + +use penf, only : I_P, R_P +use wenoof_alpha_object +use wenoof_alpha_rec_js +use wenoof_alpha_rec_m +use wenoof_alpha_rec_z +use wenoof_base_object +use wenoof_beta_object +use wenoof_beta_rec_js +use wenoof_kappa_object +use wenoof_kappa_rec_js +use wenoof_weights_object + +implicit none +private +public :: weights_js +public :: weights_js_constructor +public :: create_weights_js_constructor + +type, extends(weights_object_constructor) :: weights_js_constructor + !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights object constructor. + class(alpha_object_constructor), allocatable :: alpha_constructor !< Alpha coefficients (non linear weights) constructor. + class(beta_object_constructor), allocatable :: beta_constructor !< Beta coefficients (smoothness indicators) constructor. + class(kappa_object_constructor), allocatable :: kappa_constructor !< kappa coefficients (optimal, linear weights) constructor. +endtype weights_js_constructor + +type, extends(weights_object):: weights_js + !< Jiang-Shu and Gerolymos-Senechal-Vallet weights object. + !< + !< @note The provided WENO weights implements the weights defined in *Efficient Implementation of Weighted ENO + !< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and + !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, + !< doi:10.1016/j.jcp.2009.07.039 + class(alpha_object), allocatable :: alpha !< Alpha coefficients (non linear weights). + class(beta_object), allocatable :: beta !< Beta coefficients (smoothness indicators). + class(kappa_object), allocatable :: kappa !< kappa coefficients (optimal, linear weights). + contains + ! deferred public methods + procedure, pass(self) :: create !< Create weights. + procedure, pass(self) :: compute !< Compute weights. + procedure, pass(self) :: description !< Return weights string-description. + procedure, pass(self) :: destroy !< Destroy weights. +endtype weights_js + +contains + ! public non TBP + subroutine create_weights_js_constructor(S, alpha_constructor, beta_constructor, kappa_constructor, constructor, & + face_left, face_right, eps) + !< Create weights constructor. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(alpha_object_constructor), intent(in) :: alpha_constructor !< Alpha constructor. + class(beta_object_constructor), intent(in) :: beta_constructor !< Beta constructor. + class(kappa_object_constructor), intent(in) :: kappa_constructor !< kappa constructor. + class(weights_object_constructor), allocatable, intent(out) :: constructor !< Constructor. + logical, intent(in), optional :: face_left !< Activate left-face interpolations. + logical, intent(in), optional :: face_right !< Activate right-face interpolations. + real(R_P), intent(in), optional :: eps !< Small epsilon to avoid zero-div. + + allocate(weights_js_constructor :: constructor) + constructor%S = S + if (present(face_left)) constructor%face_left = face_left + if (present(face_right)) constructor%face_right = face_right + if (present(eps)) constructor%eps = eps + select type(constructor) + type is(weights_js_constructor) + allocate(constructor%alpha_constructor, source=alpha_constructor) + allocate(constructor%beta_constructor, source=beta_constructor) + allocate(constructor%kappa_constructor, source=kappa_constructor) + endselect + endsubroutine create_weights_js_constructor + + ! deferred public methods + subroutine create(self, constructor) + !< Create reconstructor. + class(weights_js), intent(inout) :: self !< Weights. + class(base_object_constructor), intent(in) :: constructor !< Constructor. + + call self%destroy + call self%create_(constructor=constructor) + allocate(self%values(1:2, 0:self%S - 1)) + self%values = 0._R_P + select type(constructor) + type is(weights_js_constructor) + associate(alpha_constructor=>constructor%alpha_constructor, & + beta_constructor=>constructor%beta_constructor, & + kappa_constructor=>constructor%kappa_constructor) + + select type(alpha_constructor) + type is(alpha_rec_js_constructor) + allocate(alpha_rec_js :: self%alpha) + call self%alpha%create(constructor=alpha_constructor) + type is(alpha_rec_m_constructor) + ! @TODO implement this + error stop 'alpha_rec_m to be implemented' + type is(alpha_rec_z_constructor) + ! @TODO implement this + error stop 'alpha_rec_z to be implemented' + endselect + + select type(beta_constructor) + type is(beta_rec_js_constructor) + allocate(beta_rec_js :: self%beta) + call self%beta%create(constructor=beta_constructor) + endselect + + select type(kappa_constructor) + type is(kappa_rec_js_constructor) + allocate(kappa_rec_js :: self%kappa) + call self%kappa%create(constructor=kappa_constructor) + endselect + endassociate + endselect + endsubroutine create + + pure subroutine compute(self, stencil) + !< Compute weights. + class(weights_js), intent(inout) :: self !< Weights. + real(R_P), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + integer(I_P) :: f, s !< Counters. + + call self%beta%compute(stencil=stencil) + call self%alpha%compute(beta=self%beta, kappa=self%kappa) + do s=0, self%S - 1 ! stencils loop + do f=self%f1, self%f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + self%values(f + self%ff, s) = self%alpha%values(f, s) / self%alpha%values_sum(f) + enddo + enddo + endsubroutine compute + + pure function description(self) result(string) + !< Return string-description of weights. + class(weights_js), intent(in) :: self !< Weights. + character(len=:), allocatable :: string !< String-description. + +#ifndef DEBUG + ! error stop in pure procedure is a F2015 feature not yet supported in debug mode + error stop 'weights_js%description to be implemented, do not use!' +#endif + endfunction description + + elemental subroutine destroy(self) + !< Destroy weights. + class(weights_js), intent(inout) :: self !< Weights. + + call self%destroy_ + if (allocated(self%values)) deallocate(self%values) + if (allocated(self%alpha)) deallocate(self%alpha) + if (allocated(self%beta)) deallocate(self%beta) + if (allocated(self%kappa)) deallocate(self%kappa) + endsubroutine destroy +endmodule wenoof_weights_js diff --git a/src/lib/wenoof.f90 b/src/lib/wenoof.f90 index 865f22d..174ca81 100644 --- a/src/lib/wenoof.f90 +++ b/src/lib/wenoof.f90 @@ -3,55 +3,85 @@ module wenoof !< WenOOF, WENO interpolation Object Oriented Fortran library use penf -use wenoof_alpha_coefficients -use wenoof_alpha_coefficients_js -use wenoof_alpha_coefficients_z -use wenoof_alpha_coefficients_m -use wenoof_interpolator -use wenoof_interpolator_js -use wenoof_optimal_weights -use wenoof_optimal_weights_js -use wenoof_polynomials -use wenoof_polynomials_js -use wenoof_smoothness_indicators -use wenoof_smoothness_indicators_js +use wenoof_alpha_object +use wenoof_alpha_rec_js +use wenoof_alpha_rec_z +use wenoof_alpha_rec_m +use wenoof_beta_object +use wenoof_beta_rec_js +use wenoof_interpolator_object +! use wenoof_interpolator_js +use wenoof_interpolations_object +use wenoof_interpolations_rec_js +use wenoof_kappa_object +use wenoof_kappa_rec_js +use wenoof_reconstructor_js +use wenoof_weights_object +use wenoof_weights_js implicit none private -public :: interpolator +public :: interpolator_object public :: wenoof_create contains - subroutine wenoof_create(interpolator_type, S, wenoof_interpolator, eps) + subroutine wenoof_create(interpolator_type, S, interpolator, face_left, face_right, eps) !< WenOOF creator, create a concrete WENO interpolator object. - character(*), intent(in) :: interpolator_type !< Type of the interpolator. - integer(I_P), intent(in) :: S !< Stencil dimension. - class(interpolator), allocatable, intent(out) :: wenoof_interpolator !< The concrete WENO interpolator. - real(R_P), intent(in), optional :: eps !< Parameter for avoiding division by zero. - class(smoothness_indicators_constructor), allocatable :: is_constructor !< Smoothness indicators constructor. - class(alpha_coefficients_constructor), allocatable :: alpha_constructor !< Alpha coefficients constructor. - class(optimal_weights_constructor), allocatable :: weights_constructor !< Optimal weights constructor. - class(polynomials_constructor), allocatable :: polynom_constructor !< Polynomials constructor. - class(interpolator_constructor), allocatable :: interp_constructor !< Interpolator constructor. + character(*), intent(in) :: interpolator_type !< Type of the interpolator. + integer(I_P), intent(in) :: S !< Stencil dimension. + class(interpolator_object), allocatable, intent(out) :: interpolator !< The concrete WENO interpolator. + logical, intent(in), optional :: face_left !< Activate left-face interpolations. + logical, intent(in), optional :: face_right !< Activate right-face interpolations. + real(R_P), intent(in), optional :: eps !< Small epsilon to avoid zero-div. + class(alpha_object_constructor), allocatable :: alpha_constructor !< Alpha constructor. + class(beta_object_constructor), allocatable :: beta_constructor !< Beta constructor. + class(interpolations_object_constructor), allocatable :: interpolations_constructor !< Interpolations constructor. + class(kappa_object_constructor), allocatable :: kappa_constructor !< Kappa constructor. + class(weights_object_constructor), allocatable :: weights_constructor !< Weights constructor. + class(interpolator_object_constructor), allocatable :: interpolator_constructor !< Interpolator constructor. select case(trim(adjustl(interpolator_type))) - case('JS') - call create_smoothness_indicators_js_constructor(S=S, constructor=is_constructor) - call create_alpha_coefficients_js_constructor(S=S, constructor=alpha_constructor) - call create_optimal_weights_js_constructor(S=S, constructor=weights_constructor) - call create_polynomials_js_constructor(S=S, constructor=polynom_constructor) - call create_interpolator_js_constructor(is=is_constructor, & - alpha=alpha_constructor, & - weights=weights_constructor, & - polynom=polynom_constructor, & - S=S, eps=eps, & - constructor=interp_constructor) - allocate(interpolator_js :: wenoof_interpolator) - call wenoof_interpolator%create(constructor=interp_constructor) + case('interpolator-JS') + ! @TODO implement this + error stop 'interpolator-JS to be implemented' + case('reconstructor-JS') + call create_alpha_rec_js_constructor(S=S, & + constructor=alpha_constructor, & + face_left=face_left, & + face_right=face_right, & + eps=eps) + call create_beta_rec_js_constructor(S=S, & + constructor=beta_constructor, & + face_left=face_left, & + face_right=face_right) + call create_interpolations_rec_js_constructor(S=S, & + constructor=interpolations_constructor, & + face_left=face_left, & + face_right=face_right) + call create_kappa_rec_js_constructor(S=S, constructor=kappa_constructor) + call create_weights_js_constructor(S=S, & + alpha_constructor=alpha_constructor, & + beta_constructor=beta_constructor, & + kappa_constructor=kappa_constructor, & + constructor=weights_constructor, & + face_left=face_left, & + face_right=face_right, & + eps=eps) + call create_reconstructor_js_constructor(S=S, & + interpolations_constructor=interpolations_constructor, & + weights_constructor=weights_constructor, & + constructor=interpolator_constructor, & + face_left=face_left, & + face_right=face_right, & + eps=eps) + allocate(reconstructor_js :: interpolator) + call interpolator%create(constructor=interpolator_constructor) case('JS-Z') ! @TODO add Z support + error stop 'JS-Z to be implemented' case('JS-M') ! @TODO add M support + error stop 'JS-M to be implemented' case default ! @TODO add error handling endselect diff --git a/src/lib/wenoof_alpha_coefficients.F90 b/src/lib/wenoof_alpha_coefficients.F90 deleted file mode 100644 index 387bea7..0000000 --- a/src/lib/wenoof_alpha_coefficients.F90 +++ /dev/null @@ -1,85 +0,0 @@ -!< Abstract alpha coefficients object. -module wenoof_alpha_coefficients -!< Abstract alpha coefficients object. - -use penf, only : I_P, R_P -use wenoof_base_object - -implicit none -private -public :: alpha_coefficients -public :: alpha_coefficients_constructor - -type, extends(base_object_constructor) :: alpha_coefficients_constructor - !< Abstract alpha coefficients object constructor. - integer(I_P) :: S = 0 !< Stencils dimension. -endtype alpha_coefficients_constructor - -type, extends(base_object) :: alpha_coefficients - !< Abstract alpha coefficients object. - !< - !< @note Do not implement any real alpha coefficient, but provide the interface for the different alpha coefficient implemented. - real(R_P), allocatable :: alpha_coef(:,:) !< Alpha coefficients [1:2,0:S-1] - real(R_P), allocatable :: alpha_tot(:) !< Sum of alpha coefficients - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute alpha coefficients. - procedure, nopass :: description !< Return alpha coefficients string-description. - ! public methods - procedure, pass(self) :: create !< Create alpha coefficients. - procedure, pass(self) :: destroy !< Destroy alpha coefficients. -endtype alpha_coefficients - -contains - ! deferred public methods - pure subroutine compute(self, S, weight_opt, IS, eps, f1, f2) - !< Compute alpha coefficients. - class(alpha_coefficients), intent(inout) :: self !< Alpha coefficients. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: weight_opt(1:2, 0:S-1) !< Optimal weight of the stencil. - real(R_P), intent(in) :: IS(1:2, 0:S-1) !< Smoothness indicators of the stencils. - real(R_P), intent(in) :: eps !< Parameter for avoiding divided by zero. - integer(I_P), intent(in) :: f1, f2 !< Faces to be computed. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'alpha_coefficients%compute to be implemented by your concrete alpha coefficients object' -#endif - endsubroutine compute - - pure function description() result(string) - !< Return alpha coefficients string-description. - character(len=:), allocatable :: string !< String-description. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'alpha_coefficients%description to be implemented by your concrete alpha coefficients object' -#endif - endfunction description - - ! public methods - pure subroutine create(self, constructor) - !< Create alpha coefficients. - class(alpha_coefficients), intent(inout) :: self !< Alpha coefficients. - class(base_object_constructor), intent(in) :: constructor !< Alpha coefficients constructor. - - call self%destroy - select type(constructor) - class is(alpha_coefficients_constructor) - allocate(self%alpha_coef(1:2, 0:constructor%S - 1)) - class default - ! @TODO add error handling - endselect - allocate(self%alpha_tot(1:2)) - self%alpha_coef = 0._R_P - self%alpha_tot = 0._R_P - endsubroutine create - - elemental subroutine destroy(self) - !< Destroy alpha coefficients. - class(alpha_coefficients), intent(inout) :: self !< Alpha coefficients. - - if (allocated(self%alpha_coef)) deallocate(self%alpha_coef) - if (allocated(self%alpha_tot)) deallocate(self%alpha_tot) - endsubroutine destroy -endmodule wenoof_alpha_coefficients diff --git a/src/lib/wenoof_alpha_coefficients_js.f90 b/src/lib/wenoof_alpha_coefficients_js.f90 deleted file mode 100644 index 63d5b34..0000000 --- a/src/lib/wenoof_alpha_coefficients_js.f90 +++ /dev/null @@ -1,81 +0,0 @@ -!< Jiang-Shu alpha coefficients object. -module wenoof_alpha_coefficients_js -!< Jiang-Shu alpha coefficients object. -!< -!< @note The provided WENO alpha coefficient implements the alpha coefficients defined in *Efficient Implementation of Weighted ENO -!< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. - -use penf, only : I_P, R_P -use wenoof_alpha_coefficients - -implicit none -private -public :: alpha_coefficients_js -public :: alpha_coefficients_js_constructor -public :: create_alpha_coefficients_js_constructor - -type, extends(alpha_coefficients_constructor) :: alpha_coefficients_js_constructor - !< Jiang-Shu alpha coefficient object constructor. -endtype alpha_coefficients_js_constructor - -type, extends(alpha_coefficients) :: alpha_coefficients_js - !< Jiang-Shu alpha coefficient object. - !< - !< @note The provided WENO alpha coefficient implements the alpha coefficients defined in *Efficient Implementation of Weighted - !< ENO Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. - private - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute alpha coefficients. - procedure, nopass :: description !< Return alpha coefficients string-description. -endtype alpha_coefficients_js - -contains - ! public non TBP - subroutine create_alpha_coefficients_js_constructor(S, constructor) - !< Create alpha coefficients constructor. - integer(I_P), intent(in) :: S !< Stencils dimension. - class(alpha_coefficients_constructor), allocatable, intent(out) :: constructor !< Alpha coefficients constructor. - - allocate(alpha_coefficients_js_constructor :: constructor) - constructor%S = S - endsubroutine create_alpha_coefficients_js_constructor - - ! deferred public methods - pure subroutine compute(self, S, weight_opt, IS, eps, f1, f2) - !< Compute alpha coefficients. - class(alpha_coefficients_js), intent(inout) :: self !< Alpha coefficient. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: weight_opt(1: 2, 0: S - 1) !< Optimal weight of the stencil. - real(R_P), intent(in) :: IS(1: 2, 0: S - 1) !< Smoothness indicators of the stencils. - real(R_P), intent(in) :: eps !< Parameter for avoiding divided by zero. - integer(I_P), intent(in) :: f1, f2 !< Faces to be computed. - integer(I_P) :: f, s1 !< Counters. - - self%alpha_tot = 0._R_P - do s1=0, S - 1 ! stencil loops - do f=f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - self%alpha_coef(f, s1) = weight_opt(f, s1) * (1._R_P/(eps + IS(f, s1))**S) - self%alpha_tot(f) = self%alpha_tot(f) + self%alpha_coef(f, s1) - enddo - enddo - endsubroutine compute - - pure function description() result(string) - !< Return alpha coefficients string-descripition. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - - string = 'WENO alpha coefficient'//nl - string = string//' Based on the work by Jiang and Shu "Efficient Implementation of Weighted ENO Schemes", see '// & - 'JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130'//nl - string = string//' The "compute" method has the following public API'//nl - string = string//' alpha(S, weigt_opt, IS, eps, f1, f2)'//nl - string = string//' where:'//nl - string = string//' S: integer(I_P), intent(in), the number of the stencils used'//nl - string = string//' weight_opt: real(R_P), intent(in), the optimal weight of the actual stencil'//nl - string = string//' IS: real(R_P), intent(in), the smoothness indicator of the actual stencil'//nl - string = string//' eps: real(R_P), intent(in), the coefficient to avoid zero division used'//nl - string = string//' f1, f2: integer(I_P), intent(in), the faces to be computed (1 => left interface, 2 => right interface)' - endfunction description -endmodule wenoof_alpha_coefficients_js diff --git a/src/lib/wenoof_alpha_coefficients_m.f90 b/src/lib/wenoof_alpha_coefficients_m.f90 deleted file mode 100644 index 17f917f..0000000 --- a/src/lib/wenoof_alpha_coefficients_m.f90 +++ /dev/null @@ -1,134 +0,0 @@ -!< Henrick alpha coefficients object. -module wenoof_alpha_coefficients_m -!< Henrick alpha coefficients object. -!< -!< @note The provided WENO alpha coefficient implements the alpha coefficients defined in *Mapped weighted essentially -!< non-oscillatory schemes: Achieving optimal order near critical points*, Andrew K. Henrick, Tariq D. Aslam, Joseph M. Powers, JCP, -!< 2005, vol. 207, pp. 542-567, doi:10.1016/j.jcp.2005.01.023 - -use penf, only : I_P, R_P -use wenoof_alpha_coefficients -use wenoof_alpha_coefficients_js -use wenoof_alpha_coefficients_z -use wenoof_base_object - -implicit none -private -public :: alpha_coefficients_m -public :: alpha_coefficients_m_constructor -public :: create_alpha_coefficients_m_constructor - -type, extends(alpha_coefficients_constructor) :: alpha_coefficients_m_constructor - !< Henrick alpha coefficients object constructor. - character(len=:), allocatable :: base_type !< Base alpha coefficient type. -endtype alpha_coefficients_m_constructor - -type, extends(alpha_coefficients) :: alpha_coefficients_m - !< Henrick alpha coefficients object. - !< - !< @note The provided WENO alpha coefficient implements the alpha coefficients defined in *Mapped weighted essentially - !< non-oscillatory schemes: Achieving optimal order near critical points*, Andrew K. Henrick, Tariq D. Aslam, Joseph M. Powers, - !< JCP, 2005, vol. 207, pp. 542-567, doi:10.1016/j.jcp.2005.01.023. - class(alpha_coefficients), allocatable :: alpha_base !< Base alpha coefficients to be re-mapped. - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute alpha coefficients. - procedure, nopass :: description !< Return alpha coefficients string-description. - ! public methods - procedure, pass(self) :: create !< Create alpha coefficients. - procedure, pass(self) :: destroy !< Destroy alpha coefficients. -endtype alpha_coefficients_m - -contains - ! public non TBP - subroutine create_alpha_coefficients_m_constructor(S, constructor) - !< Create alpha coefficients constructor. - !< - !< #TODO add actual M support (this is a copy of simple JS). - integer(I_P), intent(in) :: S !< Stencils dimension. - class(alpha_coefficients_constructor), allocatable, intent(out) :: constructor !< Alpha coefficients constructor. - - allocate(alpha_coefficients_m_constructor :: constructor) - constructor%S = S - endsubroutine create_alpha_coefficients_m_constructor - - ! deferred public methods - pure subroutine compute(self, S, weight_opt, IS, eps, f1, f2) - !< Compute alpha coefficients. - class(alpha_coefficients_m), intent(inout) :: self !< Alpha coefficients. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: weight_opt(1: 2, 0: S - 1) !< Optimal weight of the stencil. - real(R_P), intent(in) :: IS(1: 2, 0: S - 1) !< Smoothness indicators of the stencils. - real(R_P), intent(in) :: eps !< Parameter for avoiding divided by zero. - integer(I_P), intent(in) :: f1, f2 !< Faces to be computed. - integer(I_P) :: f, s1 !< Counters. - - self%alpha_tot = 0._R_P - call self%alpha_base%compute(S=S, weight_opt=weight_opt, IS=IS, eps=eps, f1=f1, f2=f2) - do s1=0, S - 1 ! stencil loops - do f=f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - self%alpha_coef(f, s1) = & - (self%alpha_base%alpha_coef(f, s1) * (weight_opt(f, s1) + weight_opt(f, s1) * weight_opt(f, s1) - & - 3._R_P * weight_opt(f, s1) * self%alpha_base%alpha_coef(f, s1) + self%alpha_base%alpha_coef(f, s1) * & - self%alpha_base%alpha_coef(f, s1))) / & - (weight_opt(f, s1) * weight_opt(f, s1) + self%alpha_base%alpha_coef(f, s1) * & - (1._R_P - 2._R_P * weight_opt(f, s1))) - self%alpha_tot(f) = self%alpha_tot(f) + self%alpha_coef(f, s1) - enddo - enddo - endsubroutine compute - - pure function description() result(string) - !< Return alpha coefficients string-descripition. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - - string = 'WENO alpha coefficient'//nl - string = string//' Based on the work by Henrick, Aslam and Powers "Mapped weighted essentially non-oscillatory schemes: '// & - 'Achieving optimal order near critical points", see JCP, 2005, vol. 207, pp. 542--567, '// & - 'doi:10.1006/jcph.1996.0130'//nl - string = string//' The "compute" method has the following public API'//nl - string = string//' compute(S, weigt_opt, IS, eps, f1, f2)'//nl - string = string//' where:'//nl - string = string//' S: integer(I_P), intent(in), the number of the stencils used'//nl - string = string//' weight_opt: real(R_P), intent(in), the optimal weight of the actual stencil'//nl - string = string//' IS: real(R_P), intent(in), the smoothness indicator of the actual stencil'//nl - string = string//' eps: real(R_P), intent(in), the coefficient to avoid zero division used'//nl - string = string//' f1, f2: integer(I_P), intent(in), the faces to be computed (1 => left interface, 2 => right interface)' - endfunction description - - ! overridden methods - pure subroutine create(self, constructor) - !< Create alpha coefficients. - class(alpha_coefficients_m), intent(inout) :: self !< Alpha coefficients. - class(base_object_constructor), intent(in) :: constructor !< Alpha coefficients constructor. - - call self%destroy - call self%alpha_coefficients%create(constructor=constructor) - select type(constructor) - type is(alpha_coefficients_m_constructor) - if (allocated(constructor%base_type)) then - select case(constructor%base_type) - case('JS') - if (allocated(self%alpha_base)) deallocate(self%alpha_base) - allocate(alpha_coefficients_js :: self%alpha_base) - call self%alpha_base%create(constructor=constructor) - case('Z') - if (allocated(self%alpha_base)) deallocate(self%alpha_base) - allocate(alpha_coefficients_z :: self%alpha_base) - call self%alpha_base%create(constructor=constructor) - endselect - endif - class default - ! @TODO add error handling - endselect - endsubroutine create - - elemental subroutine destroy(self) - !< Destroy alpha coefficients. - class(alpha_coefficients_m), intent(inout) :: self !< Alpha coefficients. - - call self%alpha_coefficients%destroy - if (allocated(self%alpha_base)) deallocate(self%alpha_base) - endsubroutine destroy -endmodule wenoof_alpha_coefficients_m diff --git a/src/lib/wenoof_alpha_coefficients_z.f90 b/src/lib/wenoof_alpha_coefficients_z.f90 deleted file mode 100644 index b4c9971..0000000 --- a/src/lib/wenoof_alpha_coefficients_z.f90 +++ /dev/null @@ -1,112 +0,0 @@ -!< Borges alpha coefficients object. -module wenoof_alpha_coefficients_z -!< Borges alpha coefficients object. -!< -!< @note The provided WENO alpha coefficients implements the alpha coefficients defined in *An improved weighted essentially -!< non-oscillatory scheme for hyperbolic conservation laws*, Rafael Borges, Monique Carmona, Bruno Costa and Wai Sun Don, JCP, 2008, -!< vol. 227, pp. 3191-3211, doi: 10.1016/j.jcp.2007.11.038. - -use penf, only : I_P, R_P -use wenoof_alpha_coefficients - -implicit none -private -public :: alpha_coefficients_z -public :: alpha_coefficients_z_constructor -public :: create_alpha_coefficients_z_constructor - -type, extends(alpha_coefficients_constructor) :: alpha_coefficients_z_constructor - !< Borges WENO alpha coefficients object constructor. -endtype alpha_coefficients_z_constructor - -type, extends(alpha_coefficients) :: alpha_coefficients_z - !< Borges WENO alpha coefficients object. - !< - !< @note The provided WENO alpha coefficients implements the alpha coefficients defined in *An improved weighted essentially - !< non-oscillatory scheme for hyperbolic conservation laws*, Rafael Borges, Monique Carmona, Bruno Costa and Wai Sun Don, JCP, - !< 2008, vol. 227, pp. 3191-3211, doi: 10.1016/j.jcp.2007.11.038. - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute coefficients. - procedure, nopass :: description !< Return string-description of coefficients. -endtype alpha_coefficients_z -contains - ! public non TBP - subroutine create_alpha_coefficients_z_constructor(S, constructor) - !< Create alpha coefficients constructor. - !< - !< #TODO add actual Z support (this is a copy of simple JS). - integer(I_P), intent(in) :: S !< Stencils dimension. - class(alpha_coefficients_constructor), allocatable, intent(out) :: constructor !< Alpha coefficients constructor. - - allocate(alpha_coefficients_z_constructor :: constructor) - constructor%S = S - endsubroutine create_alpha_coefficients_z_constructor - - ! deferred public methods - pure subroutine compute(self, S, weight_opt, IS, eps, f1, f2) - !< Compute alpha coefficients. - class(alpha_coefficients_z), intent(inout) :: self !< Alpha coefficients. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: weight_opt(1: 2, 0 : S - 1) !< Optimal weight of the stencil. - real(R_P), intent(in) :: IS(1: 2, 0 : S - 1) !< Smoothness indicators of the stencils. - real(R_P), intent(in) :: eps !< Parameter for avoiding divided by zero. - integer(I_P), intent(in) :: f1, f2 !< Faces to be computed. - integer(I_P) :: f, s1 !< Counters. - - self%alpha_tot = 0._R_P - do s1=0, S - 1 ! stencil loops - do f=f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - self%alpha_coef(f, s1) = weight_opt(f, s1) * ((1._R_P + (tau(S,IS)/(eps+IS(f, s1)))) ** (weno_exp(S))) - self%alpha_tot(f) = self%alpha_tot(f) + self%alpha_coef(f, s1) - enddo - enddo - endsubroutine compute - - pure function description() result(string) - !< Return alpha coefficients string-descripition. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - - string = 'WENO alpha coefficients'//nl - string = string//' Based on the work by Borges, Carmona, Costa and Don "An improved weighted essentially non-oscillatory '// & - 'scheme for hyperbolic conservation laws", see '// & - 'JCP, 2008, vol. 227, pp. 3191--3211, doi:10.1016/j.jcp.2007.11.038'//nl - string = string//' The "compute" method has the following public API'//nl - string = string//' compute(S, weigt_opt, IS, eps, f1, f2)'//nl - string = string//' where:'//nl - string = string//' S: integer(I_P), intent(in), the number of the stencils used'//nl - string = string//' weight_opt: real(R_P), intent(in), the optimal weight of the actual stencil'//nl - string = string//' IS: real(R_P), intent(in), the smoothness indicator of the actual stencil'//nl - string = string//' eps: real(R_P), intent(in), the coefficient to avoid zero division used'//nl - string = string//' f1, f2: integer(I_P), intent(in), the faces to be computed (1 => left interface, 2 => right interface)' - return - !--------------------------------------------------------------------------------------------------------------------------------- - endfunction description - - ! private non TBP - pure function tau(S, IS) result(w_tau) - !< Compute the tau coefficient used in the WENO-Z alpha coefficients. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: IS(0:S - 1) !< Smoothness indicators. - real(R_P) :: w_tau !< Tau coefficient. - - w_tau = abs(IS(0) - (1-weno_odd(S))*IS(1) - (1-weno_odd(S))*IS(S-2_I_P) + (1-2*weno_odd(S))*IS(S-1_I_P)) - endfunction tau - - pure function weno_exp(S) result(w_exp) - !< Compute the exponent used in the alpha function. - integer(I_P), intent(in) :: S !< Number of stencils used. - integer(I_P) :: w_exp !< Exponent used in the alpha function. - - w_exp = int(S, I_P) - endfunction weno_exp - - pure function weno_odd(S) result(w_odd) - !< Compute the distinguisher between odd and even number of stencils. - integer(I_P), intent(in) :: S !< Number of stencils used. - integer(I_P) :: w_odd !< Distinguishing between odd and even number of stencils. - - w_odd = int(mod(S, 2_I_P), I_P) - endfunction weno_odd -endmodule wenoof_alpha_coefficients_z diff --git a/src/lib/wenoof_base_object.f90 b/src/lib/wenoof_base_object.f90 deleted file mode 100644 index 0456b9f..0000000 --- a/src/lib/wenoof_base_object.f90 +++ /dev/null @@ -1,51 +0,0 @@ -!< Abstract base object, the ancestor of all. -module wenoof_base_object -!< Abstract base object, the ancestor of all. -!< -!< Define a minimal, base, object that is used as ancestor of all objects, e.g. smoothness indicator, optimal weights, etc... - -implicit none -private -public :: base_object -public :: base_object_constructor - -type :: base_object_constructor - !< Abstract base object constructor. -endtype base_object_constructor - -type, abstract :: base_object - !< Abstract base object, the ancestor of all. - !< - !< Define a minimal, base, object that is used as ancestor of all objects, e.g. smoothness indicator, optimal weights, etc... - contains - ! deferred public methods - procedure(create_interface), pass(self), deferred :: create !< Create object. - procedure(description_interface), nopass, deferred :: description !< Return object string-description. - procedure(destroy_interface), pass(self), deferred :: destroy !< Destroy object. -endtype base_object - -abstract interface - !< Abstract interface of [base_object] methods. - - subroutine create_interface(self, constructor) - !< Create object. - !< - !< @note Before call this method a concrete constructor must be instantiated. - import :: base_object - import :: base_object_constructor - class(base_object), intent(inout) :: self !< Object. - class(base_object_constructor), intent(in) :: constructor !< Object constructor. - endsubroutine create_interface - - pure function description_interface() result(string) - !< Return object string-description. - character(len=:), allocatable :: string !< String-description. - endfunction description_interface - - elemental subroutine destroy_interface(self) - !< Destroy object - import :: base_object - class(base_object), intent(inout) :: self !< Object. - endsubroutine destroy_interface -endinterface -endmodule wenoof_base_object diff --git a/src/lib/wenoof_interpolator.F90 b/src/lib/wenoof_interpolator.F90 deleted file mode 100644 index 11592fd..0000000 --- a/src/lib/wenoof_interpolator.F90 +++ /dev/null @@ -1,151 +0,0 @@ -!< Abstract interpolator object. -module wenoof_interpolator -!< Abstract interpolator object. - -use penf, only : I_P, R_P -use wenoof_alpha_coefficients -use wenoof_base_object -use wenoof_objects_factory -use wenoof_optimal_weights -use wenoof_smoothness_indicators -use wenoof_polynomials - -implicit none -private -public :: interpolator -public :: interpolator_constructor - -type, extends(base_object_constructor) :: interpolator_constructor - !< Abstract interpolator object constructor. - !< - !< @note Every concrete WENO interpolator implementations must define their own constructor type. - class(smoothness_indicators_constructor), allocatable :: is !< Smoothness indicators constructor. - class(alpha_coefficients_constructor), allocatable :: alpha !< Alpha coefficients constructor. - class(optimal_weights_constructor), allocatable :: weights !< Optimal weights constructor. - class(polynomials_constructor), allocatable :: polynom !< Polynomilas constructor. - contains - ! public methods - procedure, pass(self) :: create => create_interpolator_constructor !< Create interpolator constructor. - procedure, pass(self) :: destroy => destroy_interpolator_constructor !< Destroy interpolator constructor. -endtype interpolator_constructor - -type, extends(base_object) :: interpolator - !< Abstract interpolator object. - !< - !< @note Do not implement any actual interpolator: provide the interface for the different interpolators implemented. - class(smoothness_indicators), allocatable :: is !< Smoothness indicators. - class(alpha_coefficients), allocatable :: alpha !< Alpha coefficients. - class(optimal_weights), allocatable :: weights !< Optimal weights. - class(polynomials), allocatable :: polynom !< Polynomilas. - contains - ! public deferred methods - procedure, nopass :: description !< Return interpolator string-description. - procedure, pass(self) :: interpolate_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_debug !< Interpolate values (providing also debug values). - ! public methods - generic :: interpolate => interpolate_standard, interpolate_debug !< Interpolate values. - procedure, pass(self) :: create !< Create interpolator. - procedure, pass(self) :: destroy !< Destroy interpolator. -endtype interpolator - -contains - ! constructor methods - - ! public methods - subroutine create_interpolator_constructor(self, is, alpha, weights, polynom) - !< Create interpolator constructor. - class(interpolator_constructor), intent(inout) :: self !< Interpolator constructor. - class(smoothness_indicators_constructor), intent(in) :: is !< Smoothness indicators constructor. - class(alpha_coefficients_constructor), intent(in) :: alpha !< Alpha coefficients constructor. - class(optimal_weights_constructor), intent(in) :: weights !< Optimal weights constructor. - class(polynomials_constructor), intent(in) :: polynom !< Polynomilas constructor. - - call self%destroy - allocate(self%is, source=is ) - allocate(self%alpha, source=alpha ) - allocate(self%weights, source=weights) - allocate(self%polynom, source=polynom) - endsubroutine create_interpolator_constructor - - pure subroutine destroy_interpolator_constructor(self) - !< Destroy interpolator constructor. - class(interpolator_constructor), intent(inout) :: self !< Interpolator. - - if (allocated(self%is)) deallocate(self%is) - if (allocated(self%alpha)) deallocate(self%alpha) - if (allocated(self%weights)) deallocate(self%weights) - if (allocated(self%polynom)) deallocate(self%polynom) - endsubroutine destroy_interpolator_constructor - - ! interpolator methods - - ! deferred public methods - pure function description() result(string) - !< Return interpolator string-description. - character(len=:), allocatable :: string !< String-description. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'interpolator%description to be implemented by your concrete interpolator object' -#endif - endfunction description - - pure subroutine interpolate_standard(self, S, stencil, location, interpolation) - !< Interpolate values (without providing debug values). - class(interpolator), intent(inout) :: self !< Interpolator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - character(*), intent(in) :: location !< Location of interpolation: left, right, both. - real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'interpolator%interpolate to be implemented by your concrete interpolator object' -#endif - endsubroutine interpolate_standard - - pure subroutine interpolate_debug(self, S, stencil, location, interpolation, si, weights) - !< Interpolate values (providing also debug values). - class(interpolator), intent(inout) :: self !< Interpolator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - character(*), intent(in) :: location !< Location of interpolation: left, right, both. - real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. - real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'interpolator%interpolate to be implemented by your concrete interpolator object' -#endif - endsubroutine interpolate_debug - - ! public methods - subroutine create(self, constructor) - !< Create interpolator. - class(interpolator), intent(inout) :: self !< Interpolator. - class(base_object_constructor), intent(in) :: constructor !< Constructor. - type(objects_factory) :: factory !< Objects factory. - - call self%destroy - select type(constructor) - class is(interpolator_constructor) - call factory%create(constructor=constructor%is, object=self%is) - call factory%create(constructor=constructor%alpha, object=self%alpha) - call factory%create(constructor=constructor%weights, object=self%weights) - call factory%create(constructor=constructor%polynom, object=self%polynom) - class default - ! @TODO add error handling - endselect - endsubroutine create - - elemental subroutine destroy(self) - !< Destroy interpolator - class(interpolator), intent(inout) :: self !< Interpolator. - - if (allocated(self%is)) deallocate(self%is) - if (allocated(self%alpha)) deallocate(self%alpha) - if (allocated(self%weights)) deallocate(self%weights) - if (allocated(self%polynom)) deallocate(self%polynom) - endsubroutine destroy -endmodule wenoof_interpolator diff --git a/src/lib/wenoof_interpolator_js.f90 b/src/lib/wenoof_interpolator_js.f90 deleted file mode 100644 index 13c3a60..0000000 --- a/src/lib/wenoof_interpolator_js.f90 +++ /dev/null @@ -1,241 +0,0 @@ -!< Jiang-Shu (upwind) interpolator object. -module wenoof_interpolator_js -!< Jiang-Shu (upwind) interpolator object. - -use, intrinsic :: iso_fortran_env, only : stderr=>error_unit -use penf, only : I_P, R_P, str -use wenoof_alpha_coefficients -use wenoof_alpha_coefficients_m -use wenoof_alpha_coefficients_z -use wenoof_alpha_coefficients_js -use wenoof_base_object -use wenoof_interpolator -use wenoof_optimal_weights -use wenoof_optimal_weights_js -use wenoof_polynomials -use wenoof_polynomials_js -use wenoof_smoothness_indicators -use wenoof_smoothness_indicators_js - -implicit none -private -public :: interpolator_js -public :: interpolator_js_constructor -public :: create_interpolator_js_constructor - -type, extends(interpolator_constructor) :: interpolator_js_constructor - !< Jiang-Shu (upwind) interpolator object constructor. - !< - !< @note The constructed WENO interpolator implements the *Efficient Implementation of Weighted ENO Schemes*, - !< Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. - private - integer(I_P) :: S = 0 !< Stencils dimension. - real(R_P) :: eps = 10._R_P**(-6) !< Parameter for avoiding division by zero. -endtype interpolator_js_constructor - -type, extends(interpolator) :: interpolator_js - !< Jiang-Shu (upwind) interpolator object. - !< - !< @note The WENO interpolator implemented is the *Efficient Implementation of Weighted ENO Schemes*, - !< Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130. - !< - !< @note The supported accuracy formal order are: 3rd, 5th, 7th, 9th, 11th, 13th, 15th, 17th corresponding to use 2, 3, 4, 5, 6, - !< 7, 8, 9 stencils composed of 2, 3, 4, 5, 6, 7, 8, 9 values, respectively. - private - integer(I_P) :: S = 0_I_P !< Stencil dimension. - real(R_P) :: eps = 0._R_P !< Parameter for avoiding divisiion by zero. - contains - ! public deferred methods - procedure, nopass :: description !< Return interpolator string-description. - procedure, pass(self) :: interpolate_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_debug !< Interpolate values (providing also debug values). - ! public methods - procedure, pass(self) :: create !< Create interpolator. - procedure, pass(self) :: destroy !< Destroy interpolator. - ! private methods - procedure, pass(self), private :: compute_convolution !< Compute convolution. - procedure, pass(self), private :: compute_weights !< Compute weights. -endtype interpolator_js - -contains - ! public non TBP - subroutine create_interpolator_js_constructor(is, alpha, weights, polynom, S, constructor, eps) - !< Return an instance of [interpolator_js_constructor]. - class(smoothness_indicators_constructor), intent(in) :: is !< Smoothness indicators constructor. - class(alpha_coefficients_constructor), intent(in) :: alpha !< Alpha coefficients constructor. - class(optimal_weights_constructor), intent(in) :: weights !< Optimal weights constructor. - class(polynomials_constructor), intent(in) :: polynom !< Polynomilas constructor. - integer(I_P), intent(in) :: S !< Stencil dimension. - class(interpolator_constructor), allocatable, intent(out) :: constructor !< Interpolator constructor. - real(R_P), intent(in), optional :: eps !< Parameter for avoiding division by zero. - - allocate(interpolator_js_constructor :: constructor) - call constructor%create(is=is, alpha=alpha, weights=weights, polynom=polynom) - select type(constructor) - class is(interpolator_js_constructor) - constructor%S = S - if (present(eps)) constructor%eps = eps - endselect - endsubroutine create_interpolator_js_constructor - - ! public deferred methods - pure function description() result(string) - !< Return interpolator string-descripition. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - character(len=:), allocatable :: dummy_string !< Dummy string. - - string = 'Jiang-Shu WENO upwind-biased interpolator'//nl - string = string//' Based on the scheme proposed by Jiang and Shu "Efficient Implementation of Weighted ENO Schemes", see '// & - 'JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130'//nl - ! string = string//' Provide a formal order of accuracy equals to: '//trim(str(2*self%S - 1, .true.))//nl - ! string = string//' Use '//trim(str(self%S, .true.))//' stencils composed by '//trim(str(self%S, .true.))//' values'//nl - ! string = string//' The eps value used for avoiding division by zero is '//trim(str(self%eps, .true.))//nl - string = string//' The "interpolate" method has the following public API'//nl - string = string//' interpolate(S, stencil, location, interpolation)'//nl - string = string//' where:'//nl - string = string//' S: integer(I_P), intent(in), the number of stencils actually used'//nl - string = string//' stencil(1:, 1-S:-1+S): real(R_P), intent(in), the stencils used'//nl - string = string//' location: character(*), intent(in), the location of interpolation {left, right, both}'//nl - string = string//' interpolation(1:, 1-S:-1+S): realR_P, intent(out), the interpolated values'//nl - string = string//' The alpha coefficient are evaluated by the following method'//nl - ! string = string//self%alpha%description()//nl - string = string//' The smoothness indicators are evaluated by the following method'//nl - ! string = string//self%IS%description()//nl - string = string//' The polynomials are evaluated by the following method'//nl - ! string = string//self%polynom%description()//nl - string = string//' The optimal weights are evaluated by the following method'//nl - ! string = string//self%weights%description() - endfunction description - - pure subroutine interpolate_standard(self, S, stencil, location, interpolation) - !< Interpolate values (without providing debug values). - class(interpolator_js), intent(inout) :: self !< Interpolator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - character(*), intent(in) :: location !< Location of interpolation: left, right, both. - real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(R_P) :: weights(1:2, 0:S - 1) !< Weights of the stencils, [1:2, 0:S-1]. - integer(I_P) :: f1, f2, ff !< Faces to be computed. - integer(I_P) :: f, k !< Counters. - - call compute_faces_indexes(location=location, f1=f1, f2=f2, ff=ff) - - call self%polynom%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff) - - call self%is%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff) - - call self%alpha%compute(S=S, weight_opt=self%weights%opt, IS=self%IS%si, eps=self%eps, f1=f1, f2=f2) - - call self%compute_weights(S=S, f1=f1, f2=f2, ff=ff, weights=weights) - - call self%compute_convolution(S=S, f1=f1, f2=f2, ff=ff, weights=weights, interpolation=interpolation) - endsubroutine interpolate_standard - - pure subroutine interpolate_debug(self, S, stencil, location, interpolation, si, weights) - !< Interpolate values (providing also debug values). - class(interpolator_js), intent(inout) :: self !< Interpolator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - character(*), intent(in) :: location !< Location of interpolation: left, right, both. - real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. - real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - integer(I_P) :: f1, f2, ff !< Faces to be computed. - integer(I_P) :: f, k !< Counters. - - call compute_faces_indexes(location=location, f1=f1, f2=f2, ff=ff) - - call self%polynom%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff) - - call self%is%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff) - - call self%alpha%compute(S=S, weight_opt=self%weights%opt, IS=self%IS%si, eps=self%eps, f1=f1, f2=f2) - - call self%compute_weights(S=S, f1=f1, f2=f2, ff=ff, weights=weights) - - call self%compute_convolution(S=S, f1=f1, f2=f2, ff=ff, weights=weights, interpolation=interpolation) - - associate(is => self%is) - select type(is) - class is(smoothness_indicators) - do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - si(f + ff, :) = is%si(f, :) - enddo - endselect - endassociate - endsubroutine interpolate_debug - - ! overridden methods - subroutine create(self, constructor) - !< Create interpolator. - class(interpolator_js), intent(inout) :: self !< Interpolator. - class(base_object_constructor), intent(in) :: constructor !< Constructor. - - call self%destroy - call self%interpolator%create(constructor=constructor) - select type(constructor) - type is(interpolator_js_constructor) - self%S = constructor%S - self%eps = constructor%eps - endselect - endsubroutine create - - elemental subroutine destroy(self) - !< Destoy interpolator. - class(interpolator_js), intent(inout) :: self !< Interpolator. - - call self%interpolator%destroy - self%S = 0_I_P - self%eps = 0._R_P - endsubroutine destroy - - ! private methods - pure subroutine compute_convolution(self, S, f1, f2, ff, weights, interpolation) - !< Compute convolution. - class(interpolator_js), intent(in) :: self !< Interpolator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed. - real(R_P), intent(in) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - integer(I_P) :: f, k !< Counters. - - interpolation = 0._R_P - do k = 0, S - 1 ! stencils loop - do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - interpolation(f + ff) = interpolation(f + ff) + weights(f + ff, k) * self%polynom%poly(f, k) - enddo - enddo - endsubroutine compute_convolution - - pure subroutine compute_weights(self, S, f1, f2, ff, weights) - !< Compute weights. - class(interpolator_js), intent(in) :: self !< Interpolator. - integer(I_P), intent(in) :: S !< Number of stencils actually used. - integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed. - real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - integer(I_P) :: f, k !< Counters. - - do k = 0, S - 1 ! stencils loop - do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - weights(f + ff, k) = self%alpha%alpha_coef(f, k) / self%alpha%alpha_tot(f) - enddo - enddo - endsubroutine compute_weights - - ! private non TBP - pure subroutine compute_faces_indexes(location, f1, f2, ff) - !< Compute faces indexes given the queried location. - character(*), intent(in) :: location !< Location of interpolation: left, right, both. - integer(I_P), intent(out) :: f1, f2, ff !< Faces to be computed. - - select case(location) - case('both', 'b') - f1=1_I_P; f2=2_I_P; ff=0_I_P - case('left', 'l') - f1=1_I_P; f2=1_I_P; ff=0_I_P - case('right', 'r') - f1=2_I_P; f2=2_I_P; ff=-1_I_P - endselect - endsubroutine compute_faces_indexes -endmodule wenoof_interpolator_js diff --git a/src/lib/wenoof_objects_factory.f90 b/src/lib/wenoof_objects_factory.f90 index 5ba9160..723f58f 100644 --- a/src/lib/wenoof_objects_factory.f90 +++ b/src/lib/wenoof_objects_factory.f90 @@ -2,17 +2,19 @@ module wenoof_objects_factory !< Wenoof factory. -use wenoof_alpha_coefficients -use wenoof_alpha_coefficients_js -use wenoof_alpha_coefficients_m -use wenoof_alpha_coefficients_z +use wenoof_alpha_object +use wenoof_alpha_rec_js +use wenoof_alpha_rec_m +use wenoof_alpha_rec_z use wenoof_base_object -use wenoof_optimal_weights -use wenoof_optimal_weights_js -use wenoof_polynomials -use wenoof_polynomials_js -use wenoof_smoothness_indicators -use wenoof_smoothness_indicators_js +use wenoof_beta_object +use wenoof_beta_rec_js +use wenoof_kappa_object +use wenoof_kappa_rec_js +use wenoof_interpolations_object +use wenoof_interpolations_rec_js +use wenoof_weights_object +use wenoof_weights_js implicit none private @@ -22,17 +24,19 @@ module wenoof_objects_factory !< Factory, create an instance of concrete extension of [[base_object]] given its constructor. contains ! public methods - generic :: create => create_alpha_coefficients, & - create_optimal_weights, & - create_polynomials, & - create_smoothness_indicators !< Create a concrete instance of [[alpha_coefficients]] or - !< [[optimal_weights]] or [[polynomials]] or [[smoothness_indicators]]. - procedure, nopass :: create_base_object !< Create a concrete instance of [[base_object]]. + generic :: create => create_alpha_object, & + create_beta_object, & + create_kappa_object, & + create_interpolations_object, & + create_weights_object !< Create a concrete instance of [[alpha_object]], [[beta_object]], + !< [[kappa_object]], [[interpolations_object]] or [[weights_object]]. + procedure, nopass :: create_base_object !< Create a concrete instance of [[base_object]]. ! private methods - procedure, nopass, private :: create_alpha_coefficients !< Create a concrete instance of [[alpha_coefficients]]. - procedure, nopass, private :: create_optimal_weights !< Create a concrete instance of [[optimal_weights]]. - procedure, nopass, private :: create_polynomials !< Create a concrete instance of [[polynomials]]. - procedure, nopass, private :: create_smoothness_indicators !< Create a concrete instance of [[smoothness_indicators]]. + procedure, nopass, private :: create_alpha_object !< Create a concrete instance of [[alpha_object]]. + procedure, nopass, private :: create_beta_object !< Create a concrete instance of [[beta_object]]. + procedure, nopass, private :: create_kappa_object !< Create a concrete instance of [[kappa_object]]. + procedure, nopass, private :: create_interpolations_object !< Create a concrete instance of [[interpolations_object]]. + procedure, nopass, private :: create_weights_object !< Create a concrete instance of [[weights_object]]. endtype objects_factory contains @@ -42,81 +46,97 @@ subroutine create_base_object(constructor, object) class(base_object), allocatable, intent(out) :: object !< Object. select type(constructor) - type is(alpha_coefficients_js_constructor) - allocate(alpha_coefficients_js :: object) - type is(alpha_coefficients_m_constructor) - allocate(alpha_coefficients_m :: object) - type is(alpha_coefficients_z_constructor) - allocate(alpha_coefficients_z :: object) - type is(optimal_weights_js_constructor) - allocate(optimal_weights_js :: object) - type is(polynomials_js_constructor) - allocate(polynomials_js :: object) - type is(smoothness_indicators_js_constructor) - allocate(smoothness_indicators_js :: object) + type is(alpha_rec_js_constructor) + allocate(alpha_rec_js :: object) + type is(alpha_rec_m_constructor) + allocate(alpha_rec_m :: object) + type is(alpha_rec_z_constructor) + allocate(alpha_rec_z :: object) + type is(beta_rec_js_constructor) + allocate(beta_rec_js :: object) + type is(kappa_rec_js_constructor) + allocate(kappa_rec_js :: object) + type is(interpolations_rec_js_constructor) + allocate(interpolations_rec_js :: object) + type is(weights_js_constructor) + allocate(weights_js :: object) class default error stop 'error: WenOOF object factory do NOT know the constructor given' endselect call object%create(constructor=constructor) endsubroutine create_base_object - subroutine create_alpha_coefficients(constructor, object) - !< Create an instance of concrete extension of [[alpha_coefficients]] given its constructor. - class(alpha_coefficients_constructor), intent(in) :: constructor !< Constructor. - class(alpha_coefficients), allocatable, intent(out) :: object !< Object. + subroutine create_alpha_object(constructor, object) + !< Create an instance of concrete extension of [[alpha_object]] given its constructor. + class(alpha_object_constructor), intent(in) :: constructor !< Constructor. + class(alpha_object), allocatable, intent(out) :: object !< Object. select type(constructor) - type is(alpha_coefficients_js_constructor) - allocate(alpha_coefficients_js :: object) - type is(alpha_coefficients_m_constructor) - allocate(alpha_coefficients_m :: object) - type is(alpha_coefficients_z_constructor) - allocate(alpha_coefficients_z :: object) + type is(alpha_rec_js_constructor) + allocate(alpha_rec_js :: object) + type is(alpha_rec_m_constructor) + allocate(alpha_rec_m :: object) + type is(alpha_rec_z_constructor) + allocate(alpha_rec_z :: object) class default error stop 'error: WenOOF object factory do NOT know the constructor given' endselect call object%create(constructor=constructor) - endsubroutine create_alpha_coefficients + endsubroutine create_alpha_object - subroutine create_optimal_weights(constructor, object) - !< Create an instance of concrete extension of [[optimal_weights]] given its constructor. - class(optimal_weights_constructor), intent(in) :: constructor !< Constructor. - class(optimal_weights), allocatable, intent(out) :: object !< Object. + subroutine create_beta_object(constructor, object) + !< Create an instance of concrete extension of [[beta_object]] given its constructor. + class(beta_object_constructor), intent(in) :: constructor !< Constructor. + class(beta_object), allocatable, intent(out) :: object !< Object. select type(constructor) - type is(optimal_weights_js_constructor) - allocate(optimal_weights_js :: object) + type is(beta_rec_js_constructor) + allocate(beta_rec_js :: object) class default error stop 'error: WenOOF object factory do NOT know the constructor given' endselect call object%create(constructor=constructor) - endsubroutine create_optimal_weights + endsubroutine create_beta_object - subroutine create_polynomials(constructor, object) - !< Create an instance of concrete extension of [[polynomials]] given its constructor. - class(polynomials_constructor), intent(in) :: constructor !< Constructor. - class(polynomials), allocatable, intent(out) :: object !< Object. + subroutine create_kappa_object(constructor, object) + !< Create an instance of concrete extension of [[kappa_object]] given its constructor. + class(kappa_object_constructor), intent(in) :: constructor !< Constructor. + class(kappa_object), allocatable, intent(out) :: object !< Object. select type(constructor) - type is(polynomials_js_constructor) - allocate(polynomials_js :: object) + type is(kappa_rec_js_constructor) + allocate(kappa_rec_js :: object) class default error stop 'error: WenOOF object factory do NOT know the constructor given' endselect call object%create(constructor=constructor) - endsubroutine create_polynomials + endsubroutine create_kappa_object - subroutine create_smoothness_indicators(constructor, object) - !< Create an instance of concrete extension of [[smoothness_indicators]] given its constructor. - class(smoothness_indicators_constructor), intent(in) :: constructor !< Constructor. - class(smoothness_indicators), allocatable, intent(out) :: object !< Object. + subroutine create_interpolations_object(constructor, object) + !< Create an instance of concrete extension of [[interpolations_object]] given its constructor. + class(interpolations_object_constructor), intent(in) :: constructor !< Constructor. + class(interpolations_object), allocatable, intent(out) :: object !< Object. select type(constructor) - type is(smoothness_indicators_js_constructor) - allocate(smoothness_indicators_js :: object) + type is(interpolations_rec_js_constructor) + allocate(interpolations_rec_js :: object) class default error stop 'error: WenOOF object factory do NOT know the constructor given' endselect call object%create(constructor=constructor) - endsubroutine create_smoothness_indicators + endsubroutine create_interpolations_object + + subroutine create_weights_object(constructor, object) + !< Create an instance of concrete extension of [[weights_object]] given its constructor. + class(weights_object_constructor), intent(in) :: constructor !< Constructor. + class(weights_object), allocatable, intent(out) :: object !< Object. + + select type(constructor) + type is(weights_js_constructor) + allocate(weights_js :: object) + class default + error stop 'error: WenOOF object factory do NOT know the constructor given' + endselect + call object%create(constructor=constructor) + endsubroutine create_weights_object endmodule wenoof_objects_factory diff --git a/src/lib/wenoof_optimal_weights.F90 b/src/lib/wenoof_optimal_weights.F90 deleted file mode 100644 index 7a1b668..0000000 --- a/src/lib/wenoof_optimal_weights.F90 +++ /dev/null @@ -1,77 +0,0 @@ -!< Abstract optimal weights object. -module wenoof_optimal_weights -!< Abstract optimal weights object. - -use penf, only : I_P, R_P -use wenoof_base_object - -implicit none -private -public :: optimal_weights -public :: optimal_weights_constructor - -type, extends(base_object_constructor) :: optimal_weights_constructor - !< Abstract optimal weights object constructor. - integer(I_P) :: S = 0 !< Stencils dimension. -endtype optimal_weights_constructor - -type, extends(base_object) :: optimal_weights - !< Optimal weights object. - real(R_P), allocatable :: opt(:,:) !< Optimal weights [1:2,0:S-1]. - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute weights. - procedure, nopass :: description !< Return weights string-description. - ! public methods - procedure, pass(self) :: create !< Createte weights. - procedure, pass(self) :: destroy !< Destroy weights. -endtype optimal_weights - -contains - ! deferred public methods - pure subroutine compute(self, S) - !< Compute weights. - class(optimal_weights), intent(inout) :: self !< Optimal weights. - integer(I_P), intent(in) :: S !< Number of stencils used. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'optimal_weights%compute to be implemented by your concrete optimal weights object' -#endif - endsubroutine compute - - pure function description() result(string) - !< Return weights string-description. - character(len=:), allocatable :: string !< String-description. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'optimal_weights%description to be implemented by your concrete optimal weights object' -#endif - endfunction description - - ! public methods - pure subroutine create(self, constructor) - !< Create weights. - !< - !< @note During creation the weights are also computed. - class(optimal_weights), intent(inout) :: self !< Optimal weights. - class(base_object_constructor), intent(in) :: constructor !< Optimal weights constructor. - - call self%destroy - select type(constructor) - class is(optimal_weights_constructor) - allocate(self%opt(1:2, 0:constructor%S - 1)) - call self%compute(S=constructor%S) - class default - ! @TODO add error handling - endselect - endsubroutine create - - elemental subroutine destroy(self) - !< Destroy weights. - class(optimal_weights), intent(inout) :: self !< Optimial weights. - - if (allocated(self%opt)) deallocate(self%opt) - endsubroutine destroy -endmodule wenoof_optimal_weights diff --git a/src/lib/wenoof_optimal_weights_js.f90 b/src/lib/wenoof_optimal_weights_js.f90 deleted file mode 100644 index ca8a6f0..0000000 --- a/src/lib/wenoof_optimal_weights_js.f90 +++ /dev/null @@ -1,185 +0,0 @@ -!< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights. -module wenoof_optimal_weights_js -!< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights. -!< -!< @note The provided WENO optimal weights implements the optimal weights defined in *Efficient Implementation of Weighted ENO -!< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and -!< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, -!< doi:10.1016/j.jcp.2009.07.039 - -use penf, only : I_P, R_P -use wenoof_optimal_weights - -implicit none -private -public :: optimal_weights_js -public :: optimal_weights_js_constructor -public :: create_optimal_weights_js_constructor - -type, extends(optimal_weights_constructor) :: optimal_weights_js_constructor - !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights object constructor. -endtype optimal_weights_js_constructor - -type, extends(optimal_weights):: optimal_weights_js - !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights object. - !< - !< @note The provided WENO optimal weights implements the optimal weights defined in *Efficient Implementation of Weighted ENO - !< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and - !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524, - !< doi:10.1016/j.jcp.2009.07.039 - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute weights. - procedure, nopass :: description !< Return weights string-description. -endtype optimal_weights_js - -contains - ! public non TBP - subroutine create_optimal_weights_js_constructor(S, constructor) - !< Create optimal weights constructor. - integer(I_P), intent(in) :: S !< Stencils dimension. - class(optimal_weights_constructor), allocatable, intent(out) :: constructor !< Optimal weights constructor. - - allocate(optimal_weights_js_constructor :: constructor) - constructor%S = S - endsubroutine create_optimal_weights_js_constructor - - ! deferred public methods - pure subroutine compute(self, S) - !< Compute weights. - class(optimal_weights_js), intent(inout) :: self !< Optimal weights. - integer(I_P), intent(in) :: S !< Number of stencils used. - - associate(opt => self%opt) - select case(S) - case(2) ! 3rd order - ! 1 => left interface (i-1/2) - opt(1, 0) = 2._R_P/3._R_P ! stencil 0 - opt(1, 1) = 1._R_P/3._R_P ! stencil 1 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/3._R_P ! stencil 0 - opt(2, 1) = 2._R_P/3._R_P ! stencil 1 - case(3) ! 5th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 0.3_R_P ! stencil 0 - opt(1, 1) = 0.6_R_P ! stencil 1 - opt(1, 2) = 0.1_R_P ! stencil 2 - ! 2 => right interface (i+1/2) - opt(2, 0) = 0.1_R_P ! stencil 0 - opt(2, 1) = 0.6_R_P ! stencil 1 - opt(2, 2) = 0.3_R_P ! stencil 2 - case(4) ! 7th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 4._R_P/35._R_P ! stencil 0 - opt(1, 1) = 18._R_P/35._R_P ! stencil 1 - opt(1, 2) = 12._R_P/35._R_P ! stencil 2 - opt(1, 3) = 1._R_P/35._R_P ! stencil 3 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/35._R_P ! stencil 0 - opt(2, 1) = 12._R_P/35._R_P ! stencil 1 - opt(2, 2) = 18._R_P/35._R_P ! stencil 2 - opt(2, 3) = 4._R_P/35._R_P ! stencil 3 - case(5) ! 9th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 5._R_P/126._R_P ! stencil 0 - opt(1, 1) = 20._R_P/63._R_P ! stencil 1 - opt(1, 2) = 10._R_P/21._R_P ! stencil 2 - opt(1, 3) = 10._R_P/63._R_P ! stencil 3 - opt(1, 4) = 1._R_P/126._R_P ! stencil 4 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/126._R_P ! stencil 0 - opt(2, 1) = 10._R_P/63._R_P ! stencil 1 - opt(2, 2) = 10._R_P/21._R_P ! stencil 2 - opt(2, 3) = 20._R_P/63._R_P ! stencil 3 - opt(2, 4) = 5._R_P/126._R_P ! stencil 4 - case(6) ! 11th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 1._R_P/77._R_P ! stencil 0 - opt(1, 1) = 25._R_P/154._R_P ! stencil 1 - opt(1, 2) = 100._R_P/231._R_P ! stencil 2 - opt(1, 3) = 25._R_P/77._R_P ! stencil 3 - opt(1, 4) = 5._R_P/77._R_P ! stencil 4 - opt(1, 5) = 1._R_P/462._R_P ! stencil 5 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/462._R_P ! stencil 0 - opt(2, 1) = 5._R_P/77._R_P ! stencil 1 - opt(2, 2) = 25._R_P/77._R_P ! stencil 2 - opt(2, 3) = 100._R_P/231._R_P ! stencil 3 - opt(2, 4) = 25._R_P/154._R_P ! stencil 4 - opt(2, 5) = 1._R_P/77._R_P ! stencil 5 - case(7) ! 13th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 7._R_P/1716._R_P ! stencil 0 - opt(1, 1) = 21._R_P/286._R_P ! stencil 1 - opt(1, 2) = 175._R_P/572._R_P ! stencil 2 - opt(1, 3) = 175._R_P/429._R_P ! stencil 3 - opt(1, 4) = 105._R_P/572._R_P ! stencil 4 - opt(1, 5) = 7._R_P/286._R_P ! stencil 5 - opt(1, 6) = 1._R_P/1716._R_P ! stencil 6 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/1716._R_P ! stencil 0 - opt(2, 1) = 7._R_P/286._R_P ! stencil 1 - opt(2, 2) = 105._R_P/572._R_P ! stencil 2 - opt(2, 3) = 175._R_P/429._R_P ! stencil 3 - opt(2, 4) = 175._R_P/572._R_P ! stencil 4 - opt(2, 5) = 21._R_P/286._R_P ! stencil 5 - opt(2, 6) = 7._R_P/1716._R_P ! stencil 6 - case(8) ! 15th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 8._R_P/6435._R_P ! stencil 0 - opt(1, 1) = 196._R_P/6435._R_P ! stencil 1 - opt(1, 2) = 392._R_P/2145._R_P ! stencil 2 - opt(1, 3) = 490._R_P/1287._R_P ! stencil 3 - opt(1, 4) = 392._R_P/1287._R_P ! stencil 4 - opt(1, 5) = 196._R_P/2145._R_P ! stencil 5 - opt(1, 6) = 56._R_P/6435._R_P ! stencil 6 - opt(1, 7) = 1._R_P/6435._R_P ! stencil 7 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/6435._R_P ! stencil 0 - opt(2, 1) = 56._R_P/6435._R_P ! stencil 1 - opt(2, 2) = 196._R_P/2145._R_P ! stencil 2 - opt(2, 3) = 392._R_P/1287._R_P ! stencil 3 - opt(2, 4) = 490._R_P/1287._R_P ! stencil 4 - opt(2, 5) = 392._R_P/2145._R_P ! stencil 5 - opt(2, 6) = 196._R_P/6435._R_P ! stencil 6 - opt(2, 7) = 8._R_P/6435._R_P ! stencil 7 - case(9) ! 17th order - ! 1 => left interface (i-1/2) - opt(1, 0) = 9._R_P/24310._R_P ! stencil 0 - opt(1, 1) = 144._R_P/12155._R_P ! stencil 1 - opt(1, 2) = 1176._R_P/12155._R_P ! stencil 2 - opt(1, 3) = 3528._R_P/12155._R_P ! stencil 3 - opt(1, 4) = 882._R_P/2431._R_P ! stencil 4 - opt(1, 5) = 2352._R_P/12155._R_P ! stencil 5 - opt(1, 6) = 504._R_P/12155._R_P ! stencil 6 - opt(1, 7) = 36._R_P/12155._R_P ! stencil 7 - opt(1, 8) = 1._R_P/24310._R_P ! stencil 8 - ! 2 => right interface (i+1/2) - opt(2, 0) = 1._R_P/24310._R_P ! stencil 0 - opt(2, 1) = 36._R_P/12155._R_P ! stencil 1 - opt(2, 2) = 504._R_P/12155._R_P ! stencil 2 - opt(2, 3) = 2352._R_P/12155._R_P ! stencil 3 - opt(2, 4) = 882._R_P/2431._R_P ! stencil 4 - opt(2, 5) = 3528._R_P/12155._R_P ! stencil 5 - opt(2, 6) = 1176._R_P/12155._R_P ! stencil 6 - opt(2, 7) = 144._R_P/12155._R_P ! stencil 7 - opt(2, 8) = 9._R_P/24310._R_P ! stencil 8 - endselect - endassociate - endsubroutine compute - - pure function description() result(string) - !< Return string-description of weights. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line character. - - string = 'WENO optimal weights,'//nl - string = string//' based on the work by Jiang and Shu "Efficient Implementation of Weighted ENO Schemes", see '// & - 'JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and'//nl - string = string//' on the work by Gerolymos, Senechal and Vallet "Very-high-order weno schemes", see '// & - 'JCP, 2009, vol. 228, pp. 8481--8524, doi:10.1016/j.jcp.2009.07.039'//nl - string = string//' The optimal weights are allocated in a two-dimensional array, in which the first index'//nl - string = string//' is the face selected (1 => i-1/2, 2 => i+1/2) and the second index is the number of the stencil '//nl - string = string//' (from 0 to S-1)' - endfunction description -endmodule wenoof_optimal_weights_js diff --git a/src/lib/wenoof_polynomials.F90 b/src/lib/wenoof_polynomials.F90 deleted file mode 100644 index 5ac8c31..0000000 --- a/src/lib/wenoof_polynomials.F90 +++ /dev/null @@ -1,77 +0,0 @@ -!< Abstract polynomials object. -module wenoof_polynomials -!< Abstract polynomials object. - -use penf, only : I_P, R_P -use wenoof_base_object - -implicit none -private -public :: polynomials -public :: polynomials_constructor - -type, extends(base_object_constructor) :: polynomials_constructor - !< Abstract polynomials object constructor. - integer(I_P) :: S = 0 !< Stencils dimension. -endtype polynomials_constructor - -type, extends(base_object) :: polynomials - !< Abstract polynomials object. - real(R_P), allocatable :: poly(:,:) !< Polynomial reconstructions [1:2,0:S-1]. - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute polynomials. - procedure, nopass :: description !< Return polynomials string-description. - ! public methods - procedure, pass(self) :: create !< Createte polynomials. - procedure, pass(self) :: destroy !< Destroy polynomials. -endtype polynomials - -contains - ! deferred public methods - pure subroutine compute(self, S, stencil, f1, f2, ff) - !< Compute polynomials. - class(polynomials), intent(inout) :: self !< Polynomials. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'polynomials%compute to be implemented by your concrete polynomials object' -#endif - endsubroutine compute - - pure function description() result(string) - !< Return polynomials string-description. - character(len=:), allocatable :: string !< String-description. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'polynomials%description to be implemented by your concrete polynomials object' -#endif - endfunction description - - ! public methods - pure subroutine create(self, constructor) - !< Create polynomials. - class(polynomials), intent(inout) :: self !< Polynomials. - class(base_object_constructor), intent(in) :: constructor !< Polynomials constructor. - - call self%destroy - select type(constructor) - class is(polynomials_constructor) - allocate(self%poly(1:2, 0:constructor%S - 1)) - class default - ! @TODO add error handling - endselect - self%poly = 0._R_P - endsubroutine create - - elemental subroutine destroy(self) - !< Destroy polynomials. - class(polynomials), intent(inout) :: self !< Polynomials. - - if (allocated(self%poly)) deallocate(self%poly) - endsubroutine destroy -endmodule wenoof_polynomials diff --git a/src/lib/wenoof_smoothness_indicators.F90 b/src/lib/wenoof_smoothness_indicators.F90 deleted file mode 100644 index c76d90e..0000000 --- a/src/lib/wenoof_smoothness_indicators.F90 +++ /dev/null @@ -1,77 +0,0 @@ -!< Abstract smoothness indicator object. -module wenoof_smoothness_indicators -!< Abstract smoothness indicator object. - -use penf, only : I_P, R_P -use wenoof_base_object - -implicit none -private -public :: smoothness_indicators -public :: smoothness_indicators_constructor - -type, extends(base_object_constructor) :: smoothness_indicators_constructor - !< Abstract smoothness indicators object constructor. - integer(I_P) :: S = 0 !< Stencils dimension. -endtype smoothness_indicators_constructor - -type, extends(base_object) :: smoothness_indicators - !< Abstract smoothness indicator object. - real(R_P), allocatable :: si(:,:) !< Smoothness indicators [1:2,0:S-1]. - contains - ! deferred public methods - procedure, pass(self) :: compute !< Compute smoothness indicators. - procedure, nopass :: description !< Return smoothness indicators string-description. - ! public methods - procedure, pass(self) :: create !< Createte smoothness indicators. - procedure, pass(self) :: destroy !< Destroy smoothness indicators. -endtype smoothness_indicators - -contains - ! deferred public methods - pure subroutine compute(self, S, stencil, f1, f2, ff) - !< Compute smoothness indicators. - class(smoothness_indicators), intent(inout) :: self !< Smoothness indicators. - integer(I_P), intent(in) :: S !< Number of stencils used. - real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'smoothness_indicators%compute to be implemented by your concrete smoothness indicators object' -#endif - endsubroutine compute - - pure function description() result(string) - !< Return smoothness indicators string-description. - character(len=:), allocatable :: string !< String-description. - -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'smoothness_indicators%description to be implemented by your concrete smoothness indicators object' -#endif - endfunction description - - ! public methods - pure subroutine create(self, constructor) - !< Create smoothness indicators. - class(smoothness_indicators), intent(inout) :: self !< Smoothness indicators. - class(base_object_constructor), intent(in) :: constructor !< Smoothness indicators constructor. - - call self%destroy - select type(constructor) - class is(smoothness_indicators_constructor) - allocate(self%si(1:2, 0:constructor%S - 1)) - class default - ! @TODO add error handling - endselect - self%si = 0._R_P - endsubroutine create - - elemental subroutine destroy(self) - !< Destroy smoothness indicators. - class(smoothness_indicators), intent(inout) :: self !< Smoothness indicators. - - if (allocated(self%si)) deallocate(self%si) - endsubroutine destroy -endmodule wenoof_smoothness_indicators diff --git a/src/tests/sin_reconstruction.f90 b/src/tests/sin_reconstruction.f90 index 6732d75..d34d6e5 100644 --- a/src/tests/sin_reconstruction.f90 +++ b/src/tests/sin_reconstruction.f90 @@ -6,16 +6,16 @@ module test_module use flap, only : command_line_interface use penf, only : I_P, R_P, FR_P, str, strz use pyplot_module, only : pyplot -use wenoof, only : interpolator, wenoof_create +use wenoof, only : interpolator_object, wenoof_create implicit none private public :: test -character(99), parameter :: interpolators(1:4) = ["all ", & - "JS ", & - "JS-Z", & - "JS-M"] !< List of available interpolators. +character(99), parameter :: interpolators(1:4) = ["all ", & + "reconstructor-JS", & + "JS-Z ", & + "JS-M "] !< List of available interpolators. real(R_P), parameter :: pi = 4._R_P * atan(1._R_P) !< Pi greek. type :: solution_data @@ -138,7 +138,7 @@ subroutine compute_reference_solution(self) self%solution(pn, s)%Dx = 2 * pi / self%points_number(pn) ! compute the values used for the interpolation/reconstruction of sin function: cell values do i=1 - self%S(s), self%points_number(pn) + self%S(s) - self%solution(pn, s)%x_cell(i) = i * self%solution(pn, s)%Dx + self%solution(pn, s)%x_cell(i) = i * self%solution(pn, s)%Dx - self%solution(pn, s)%Dx / 2._R_P self%solution(pn, s)%fx_cell(i) = sin(self%solution(pn, s)%x_cell(i)) enddo ! values to which the interpolation/reconstruction should tend @@ -180,7 +180,8 @@ subroutine set_cli() "sin_reconstruction --interpolator JS-Z -r ", & "sin_reconstruction --interpolator JS-M ", & "sin_reconstruction --interpolator all -p -r "]) - call cli%add(switch='--interpolator', switch_ab='-i', help='WENO interpolator type', required=.false., def='JS', act='store') + call cli%add(switch='--interpolator', switch_ab='-i', help='WENO interpolator type', required=.false., & + def='reconstructor-JS', act='store', choices='all,reconstructor-JS') call cli%add(switch='--points_number', switch_ab='-pn', nargs='+', help='Number of points used to discretize the domain', & required=.false., act='store', def='50') call cli%add(switch='--stencils', switch_ab='-s', nargs='+', help='Stencils dimensions (and number)', & @@ -208,66 +209,35 @@ subroutine parse_cli() call self%cli%get(switch='-p', val=self%plots, error=self%error) ; if (self%error/=0) stop call self%cli%get(switch='--output', val=self%output_bname, error=self%error) ; if (self%error/=0) stop call self%cli%get(switch='--errors_analysis', val=self%errors_analysis, error=self%error) ; if (self%error/=0) stop - - if (.not.is_interpolator_valid()) then - print "(A)", 'error: the interpolator type "'//trim(adjustl(self%interpolator_type))//'" is unknown!' - print "(A)", list_interpolators() - stop - endif endsubroutine parse_cli - - function is_interpolator_valid() - !< Verify if the selected interpolator is valid. - logical :: is_interpolator_valid !< Return true is the selected interpolator is available. - integer(I_P) :: s !< Counter. - - is_interpolator_valid = .false. - do s=1, size(interpolators, dim=1) - is_interpolator_valid = (trim(adjustl(self%interpolator_type))==trim(adjustl(interpolators(s)))) - if (is_interpolator_valid) exit - enddo - endfunction is_interpolator_valid - - function list_interpolators() result(list) - !< List available solvers. - character(len=:), allocatable :: list !< Pretty printed list of available interpolators. - integer(I_P) :: s !< Counter. - - list = 'Valid interpolator names are:' // new_line('a') - do s=1, ubound(interpolators, dim=1) - list = list // ' + ' // trim(adjustl(interpolators(s))) // new_line('a') - enddo - endfunction list_interpolators endsubroutine initialize subroutine perform(self) !< Perform the test. - class(test), intent(inout) :: self !< Test. - real(R_P), allocatable :: error(:,:) !< Error (norm L2) with respect the exact solution. - real(R_P), allocatable :: order(:,:) !< Observed order based on subsequent refined solutions. - class(interpolator), allocatable :: weno_interpolator !< WENO interpolator. - real(R_P), allocatable :: stencil(:,:) !< Stencils used. - integer(I_P) :: s !< Counter. - integer(I_P) :: pn !< Counter. - integer(I_P) :: i !< Counter. + class(test), intent(inout) :: self !< Test. + real(R_P), allocatable :: error(:,:) !< Error (norm L2) with respect the exact solution. + real(R_P), allocatable :: order(:,:) !< Observed order based on subsequent refined solutions. + class(interpolator_object), allocatable :: interpolator !< WENO interpolator. + real(R_P), allocatable :: stencil(:,:) !< Stencils used. + integer(I_P) :: s !< Counter. + integer(I_P) :: pn !< Counter. + integer(I_P) :: i !< Counter. call self%compute_reference_solution do s=1, self%S_number call wenoof_create(interpolator_type=trim(adjustl(self%interpolator_type)), & S=self%S(s), & - eps=self%eps, & - wenoof_interpolator=weno_interpolator) + interpolator=interpolator, & + eps=self%eps) allocate(stencil(1:2, 1-self%S(s):-1+self%S(s))) do pn=1, self%pn_number do i=1, self%points_number(pn) - stencil(1,:) = self%solution(pn, s)%fx_cell(i+0-self%S(s):i-2+self%S(s)) + stencil(1,:) = self%solution(pn, s)%fx_cell(i+1-self%S(s):i-1+self%S(s)) stencil(2,:) = self%solution(pn, s)%fx_cell(i+1-self%S(s):i-1+self%S(s)) - call weno_interpolator%interpolate(S=self%S(s), & - stencil=stencil, & - location='b', & - interpolation=self%solution(pn, s)%reconstruction(:,i), & - si=self%solution(pn, s)%si(:, i, 0:self%S(s)-1), & - weights=self%solution(pn, s)%weights(:, i, 0:self%S(s)-1)) + call interpolator%interpolate(stencil=stencil, & + interpolation=self%solution(pn, s)%reconstruction(:,i), & + si=self%solution(pn, s)%si(:, i, 0:self%S(s)-1), & + weights=self%solution(pn, s)%weights(:, i, 0:self%S(s)-1)) enddo enddo deallocate(stencil)