From 166699147253a2f72bbbb94da8275d1b42bc3e44 Mon Sep 17 00:00:00 2001 From: Stefano Zaghi Date: Fri, 7 Apr 2017 17:01:09 +0200 Subject: [PATCH 1/3] amend fobos --- fobos | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/fobos b/fobos index d84e270..f16e0c3 100644 --- a/fobos +++ b/fobos @@ -26,28 +26,28 @@ $EXDIRS = FACE/src/tests FACE/src/third_party/ # GNU [shared-gnu] template = template-shared-gnu -target = wenoof.f90 +target = wenoof.F90 build_dir = ./shared/ output = libwenoof.so mklib = shared [static-gnu] template = template-static-gnu -target = wenoof.f90 +target = wenoof.F90 build_dir = ./static/ output = libwenoof.a mklib = static [shared-gnu-debug] template = template-shared-gnu-debug -target = wenoof.f90 +target = wenoof.F90 build_dir = ./shared/ output = libwenoof.so mklib = shared [static-gnu-debug] template = template-static-gnu-debug -target = wenoof.f90 +target = wenoof.F90 build_dir = ./static/ output = libwenoof.a mklib = static @@ -63,28 +63,28 @@ build_dir = ./exe/ # Intel [shared-intel] template = template-shared-intel -target = wenoof.f90 +target = wenoof.F90 build_dir = ./shared/ output = libwenoof.so mklib = shared [static-intel] template = template-static-intel -target = wenoof.f90 +target = wenoof.F90 build_dir = ./static/ output = libwenoof.a mklib = static [shared-intel-debug] template = template-shared-intel-debug -target = wenoof.f90 +target = wenoof.F90 build_dir = ./shared/ output = libwenoof.so mklib = shared [static-intel-debug] template = template-static-intel-debug -target = wenoof.f90 +target = wenoof.F90 build_dir = ./static/ output = libwenoof.a mklib = static From b2ce4ac07ca8838473cac06968d0ac9bd615fa15 Mon Sep 17 00:00:00 2001 From: Stefano Zaghi Date: Fri, 7 Apr 2017 18:21:28 +0200 Subject: [PATCH 2/3] refactor start, do not compile --- .../wenoof_alpha_int_m.F90 | 0 .../wenoof_alpha_int_z.F90 | 0 .../wenoof_alpha_rec_m.F90 | 0 .../wenoof_alpha_rec_z.F90 | 0 .../abstract_objects/wenoof_alpha_object.F90 | 33 ++++++---- .../abstract_objects/wenoof_beta_object.F90 | 21 ++++--- .../wenoof_interpolations_object.F90 | 22 +++---- .../wenoof_interpolator_object.F90 | 32 +++++----- .../wenoof_weights_object.F90 | 23 ++++--- .../concrete_objects/wenoof_alpha_int_js.F90 | 42 +++++++------ .../concrete_objects/wenoof_alpha_rec_js.F90 | 60 +++++++++---------- .../concrete_objects/wenoof_beta_int_js.F90 | 23 +++---- .../concrete_objects/wenoof_beta_rec_js.F90 | 50 +++++++--------- .../wenoof_interpolations_int_js.F90 | 23 +++---- .../wenoof_interpolations_rec_js.F90 | 38 ++++++------ .../wenoof_interpolator_js.F90 | 30 +++++----- .../concrete_objects/wenoof_kappa_int_js.F90 | 2 +- .../wenoof_reconstructor_js.F90 | 51 ++++++++-------- .../wenoof_weights_int_js.F90 | 23 +++---- .../wenoof_weights_rec_js.F90 | 37 ++++++------ src/lib/wenoof.F90 | 6 +- 21 files changed, 263 insertions(+), 253 deletions(-) rename {src/lib/concrete_objects => quarantena}/wenoof_alpha_int_m.F90 (100%) rename {src/lib/concrete_objects => quarantena}/wenoof_alpha_int_z.F90 (100%) rename {src/lib/concrete_objects => quarantena}/wenoof_alpha_rec_m.F90 (100%) rename {src/lib/concrete_objects => quarantena}/wenoof_alpha_rec_z.F90 (100%) diff --git a/src/lib/concrete_objects/wenoof_alpha_int_m.F90 b/quarantena/wenoof_alpha_int_m.F90 similarity index 100% rename from src/lib/concrete_objects/wenoof_alpha_int_m.F90 rename to quarantena/wenoof_alpha_int_m.F90 diff --git a/src/lib/concrete_objects/wenoof_alpha_int_z.F90 b/quarantena/wenoof_alpha_int_z.F90 similarity index 100% rename from src/lib/concrete_objects/wenoof_alpha_int_z.F90 rename to quarantena/wenoof_alpha_int_z.F90 diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_m.F90 b/quarantena/wenoof_alpha_rec_m.F90 similarity index 100% rename from src/lib/concrete_objects/wenoof_alpha_rec_m.F90 rename to quarantena/wenoof_alpha_rec_m.F90 diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_z.F90 b/quarantena/wenoof_alpha_rec_z.F90 similarity index 100% rename from src/lib/concrete_objects/wenoof_alpha_rec_z.F90 rename to quarantena/wenoof_alpha_rec_z.F90 diff --git a/src/lib/abstract_objects/wenoof_alpha_object.F90 b/src/lib/abstract_objects/wenoof_alpha_object.F90 index ba3613f..da24594 100644 --- a/src/lib/abstract_objects/wenoof_alpha_object.F90 +++ b/src/lib/abstract_objects/wenoof_alpha_object.F90 @@ -23,24 +23,35 @@ module wenoof_alpha_object type, extends(base_object), abstract :: alpha_object !< Abstract alpha (non linear weights) object. - real(RPP), allocatable :: values_rank_1(:) !< Alpha values [0:S-1]. - real(RPP) :: values_sum_rank_1 !< Sum of alpha coefficients. - real(RPP), allocatable :: values_rank_2(:,:) !< Alpha values [1:2,0:S-1]. - real(RPP), allocatable :: values_sum_rank_2(:) !< Sum of alpha coefficients [1:2]. + ! real(RPP), allocatable :: values_rank_1(:) !< Alpha values [0:S-1]. + ! real(RPP) :: values_sum_rank_1 !< Sum of alpha coefficients. + ! real(RPP), allocatable :: values_rank_2(:,:) !< Alpha values [1:2,0:S-1]. + ! real(RPP), allocatable :: values_sum_rank_2(:) !< Sum of alpha coefficients [1:2]. contains ! public deferred methods - procedure(compute_interface), pass(self), deferred :: compute !< Compute alpha. + procedure(compute_interpolate_interface), pass(self), deferred :: compute_interpolate !< Compute alpha (interpolate). + procedure(compute_reconstruct_interface), pass(self), deferred :: compute_reconstruct !< Compute alpha (reconstruct). + generic :: compute => compute_interpolate, compute_reconstruct endtype alpha_object abstract interface !< Abstract interfaces of [[alpha_object]]. - pure subroutine compute_interface(self, beta, kappa) + pure subroutine compute_interpolate_interface(self, beta, kappa, values) !< 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 + class(alpha_object), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + endsubroutine compute_interpolate_interface + pure subroutine compute_reconstruct_interface(self, beta, kappa, values) + !< Compute alpha. + import :: alpha_object, beta_object, kappa_object + class(alpha_object), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + endsubroutine compute_reconstruct_interface +endinterface endmodule wenoof_alpha_object diff --git a/src/lib/abstract_objects/wenoof_beta_object.F90 b/src/lib/abstract_objects/wenoof_beta_object.F90 index 768207b..f15a9ad 100644 --- a/src/lib/abstract_objects/wenoof_beta_object.F90 +++ b/src/lib/abstract_objects/wenoof_beta_object.F90 @@ -21,30 +21,29 @@ module wenoof_beta_object type, extends(base_object), abstract :: beta_object !< Abstract Beta coefficients (smoothness indicators of stencil interpolations) object. real(RPP), allocatable :: values_rank_1(:) !< Beta values [0:S-1]. - real(RPP), allocatable :: values_rank_2(:,:) !< Beta values [1:2,0:S-1]. contains ! public methods - generic :: compute => compute_with_stencil_of_rank_1, compute_with_stencil_of_rank_2 + generic :: compute => compute_stencil_rank_1, compute_stencil_rank_2 ! deferred public methods - procedure(compute_with_stencil_of_rank_1_interface), pass(self), deferred :: compute_with_stencil_of_rank_1!< Compute beta. - procedure(compute_with_stencil_of_rank_2_interface), pass(self), deferred :: compute_with_stencil_of_rank_2!< Compute beta. + procedure(compute_stencil_rank_1_interface), pass(self), deferred :: compute_stencil_rank_1 !< Compute beta. + procedure(compute_stencil_rank_2_interface), pass(self), deferred :: compute_stencil_rank_2 !< Compute beta. endtype beta_object abstract interface !< Abstract interfaces of [[beta_object]]. - pure subroutine compute_with_stencil_of_rank_1_interface(self, stencil) + pure subroutine compute_stencil_rank_1_interface(self, stencil) !< Compute beta. import :: beta_object, RPP class(beta_object), intent(inout) :: self !< Beta. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - endsubroutine compute_with_stencil_of_rank_1_interface + endsubroutine compute_stencil_rank_1_interface - pure subroutine compute_with_stencil_of_rank_2_interface(self, stencil) + pure subroutine compute_stencil_rank_2_interface(self, stencil, values) !< Compute beta. import :: beta_object, RPP - class(beta_object), intent(inout) :: self !< Beta. - real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - endsubroutine compute_with_stencil_of_rank_2_interface + class(beta_object), intent(in) :: self !< Beta. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:,0:) !< Beta values [1:2,0:S-1]. + endsubroutine compute_stencil_rank_2_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 index 998666e..cbe349f 100644 --- a/src/lib/abstract_objects/wenoof_interpolations_object.F90 +++ b/src/lib/abstract_objects/wenoof_interpolations_object.F90 @@ -23,30 +23,30 @@ module wenoof_interpolations_object type, extends(base_object), abstract :: interpolations_object !< Abstract interpolations object. real(RPP), allocatable :: values_rank_1(:) !< Stencil interpolations values [0:S-1]. - real(RPP), allocatable :: values_rank_2(:,:) !< Stencil interpolations values [1:2,0:S-1]. + ! real(RPP), allocatable :: values_rank_2(:,:) !< Stencil interpolations values [1:2,0:S-1]. contains ! public methods - generic :: compute => compute_with_stencil_of_rank_1, compute_with_stencil_of_rank_2 + generic :: compute => compute_stencil_rank_1, compute_stencil_rank_2 ! deferred public methods - procedure(compute_with_stencil_of_rank_1_interface), pass(self), deferred :: compute_with_stencil_of_rank_1!< Compute interp. - procedure(compute_with_stencil_of_rank_2_interface), pass(self), deferred :: compute_with_stencil_of_rank_2!< Compute interp. + procedure(compute_stencil_rank_1_interface), pass(self), deferred :: compute_stencil_rank_1 !< Compute interpolation values. + procedure(compute_stencil_rank_2_interface), pass(self), deferred :: compute_stencil_rank_2 !< Compute interpolation values. endtype interpolations_object abstract interface !< Abstract interfaces of [[interpolations_object]]. - pure subroutine compute_with_stencil_of_rank_1_interface(self, stencil) + pure subroutine compute_stencil_rank_1_interface(self, stencil) !< Compute interpolations. import :: interpolations_object, RPP class(interpolations_object), intent(inout) :: self !< Interpolations. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - endsubroutine compute_with_stencil_of_rank_1_interface + endsubroutine compute_stencil_rank_1_interface - pure subroutine compute_with_stencil_of_rank_2_interface(self, stencil) + pure subroutine compute_stencil_rank_2_interface(self, stencil, values) !< Compute interpolations. import :: interpolations_object, RPP - class(interpolations_object), intent(inout) :: self !< Interpolations. - real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - endsubroutine compute_with_stencil_of_rank_2_interface + class(interpolations_object), intent(in) :: self !< Interpolations. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:, 0:) !< Interpolations values. + endsubroutine compute_stencil_rank_2_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 index 7d05c4e..d40cbc7 100644 --- a/src/lib/abstract_objects/wenoof_interpolator_object.F90 +++ b/src/lib/abstract_objects/wenoof_interpolator_object.F90 @@ -32,22 +32,18 @@ module wenoof_interpolator_object class(weights_object), allocatable :: weights !< Weights of interpolations. contains ! public methods - generic :: interpolate => interpolate_with_stencil_of_rank_1_debug, interpolate_with_stencil_of_rank_2_debug, & - interpolate_with_stencil_of_rank_1_standard, interpolate_with_stencil_of_rank_2_standard + generic :: interpolate => interpolate_stencil_rank_1_debug, interpolate_stencil_rank_2_debug, & + interpolate_stencil_rank_1_standard, interpolate_stencil_rank_2_standard ! public deferred methods - procedure(interpolate_with_stencil_of_rank_1_debug_interface), pass(self), & - deferred :: interpolate_with_stencil_of_rank_1_debug - procedure(interpolate_with_stencil_of_rank_2_debug_interface), pass(self), & - deferred :: interpolate_with_stencil_of_rank_2_debug - procedure(interpolate_with_stencil_of_rank_1_standard_interface), pass(self), & - deferred :: interpolate_with_stencil_of_rank_1_standard - procedure(interpolate_with_stencil_of_rank_2_standard_interface), pass(self), & - deferred :: interpolate_with_stencil_of_rank_2_standard + procedure(interpolate_stencil_rank_1_debug_interface), pass(self), deferred :: interpolate_stencil_rank_1_debug + procedure(interpolate_stencil_rank_2_debug_interface), pass(self), deferred :: interpolate_stencil_rank_2_debug + procedure(interpolate_stencil_rank_1_standard_interface), pass(self), deferred :: interpolate_stencil_rank_1_standard + procedure(interpolate_stencil_rank_2_standard_interface), pass(self), deferred :: interpolate_stencil_rank_2_standard endtype interpolator_object abstract interface !< Abstract interfaces of [[interpolator_object]]. - pure subroutine interpolate_with_stencil_of_rank_1_debug_interface(self, stencil, interpolation, si, weights) + pure subroutine interpolate_stencil_rank_1_debug_interface(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). import :: interpolator_object, RPP class(interpolator_object), intent(inout) :: self !< Interpolator. @@ -55,9 +51,9 @@ pure subroutine interpolate_with_stencil_of_rank_1_debug_interface(self, stencil real(RPP), intent(out) :: interpolation !< Result of the interpolation. real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [0:S-1]. real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [0:S-1]. - endsubroutine interpolate_with_stencil_of_rank_1_debug_interface + endsubroutine interpolate_stencil_rank_1_debug_interface - pure subroutine interpolate_with_stencil_of_rank_2_debug_interface(self, stencil, interpolation, si, weights) + pure subroutine interpolate_stencil_rank_2_debug_interface(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). import :: interpolator_object, RPP class(interpolator_object), intent(inout) :: self !< Interpolator. @@ -65,23 +61,23 @@ pure subroutine interpolate_with_stencil_of_rank_2_debug_interface(self, stencil real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - endsubroutine interpolate_with_stencil_of_rank_2_debug_interface + endsubroutine interpolate_stencil_rank_2_debug_interface - pure subroutine interpolate_with_stencil_of_rank_1_standard_interface(self, stencil, interpolation) + pure subroutine interpolate_stencil_rank_1_standard_interface(self, stencil, interpolation) !< Interpolate values (without providing debug values). import :: interpolator_object, RPP class(interpolator_object), intent(inout) :: self !< Interpolator. real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. real(RPP), intent(out) :: interpolation !< Result of the interpolation. - endsubroutine interpolate_with_stencil_of_rank_1_standard_interface + endsubroutine interpolate_stencil_rank_1_standard_interface - pure subroutine interpolate_with_stencil_of_rank_2_standard_interface(self, stencil, interpolation) + pure subroutine interpolate_stencil_rank_2_standard_interface(self, stencil, interpolation) !< Interpolate values (without providing debug values). import :: interpolator_object, RPP class(interpolator_object), intent(inout) :: self !< Interpolator. real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - endsubroutine interpolate_with_stencil_of_rank_2_standard_interface + endsubroutine interpolate_stencil_rank_2_standard_interface endinterface endmodule wenoof_interpolator_object diff --git a/src/lib/abstract_objects/wenoof_weights_object.F90 b/src/lib/abstract_objects/wenoof_weights_object.F90 index 9d4c8be..973656a 100644 --- a/src/lib/abstract_objects/wenoof_weights_object.F90 +++ b/src/lib/abstract_objects/wenoof_weights_object.F90 @@ -21,33 +21,33 @@ module wenoof_weights_object type, extends(base_object), abstract :: weights_object !< Weights of stencil interpolations object. real(RPP), allocatable :: values_rank_1(:) !< Weights values of stencil interpolations [0:S-1]. - real(RPP), allocatable :: values_rank_2(:,:) !< Weights values of stencil interpolations [1:2,0:S-1]. contains ! public methods - generic :: compute => compute_with_stencil_of_rank_1, compute_with_stencil_of_rank_2 - generic :: smoothness_indicators => smoothness_indicators_of_rank_1, smoothness_indicators_of_rank_2 + generic :: compute => compute_stencil_rank_1, compute_stencil_rank_2 !< Compute weights. + generic :: smoothness_indicators => smoothness_indicators_of_rank_1, smoothness_indicators_of_rank_2 !< Return IS. ! deferred public methods - procedure(compute_with_stencil_of_rank_1_interface), pass(self), deferred :: compute_with_stencil_of_rank_1 !< Compute beta. - procedure(compute_with_stencil_of_rank_2_interface), pass(self), deferred :: compute_with_stencil_of_rank_2 !< Compute beta. + procedure(compute_stencil_rank_1_interface), pass(self), deferred :: compute_stencil_rank_1 !< Compute weights. + procedure(compute_stencil_rank_2_interface), pass(self), deferred :: compute_stencil_rank_2 !< Compute weights. procedure(smoothness_indicators_of_rank_1_interface), pass(self), deferred :: smoothness_indicators_of_rank_1 !< Return IS. procedure(smoothness_indicators_of_rank_2_interface), pass(self), deferred :: smoothness_indicators_of_rank_2 !< Return IS. endtype weights_object abstract interface !< Abstract interfaces of [[weights_object]]. - pure subroutine compute_with_stencil_of_rank_1_interface(self, stencil) + pure subroutine compute_stencil_rank_1_interface(self, stencil) !< Compute beta. import :: weights_object, RPP class(weights_object), intent(inout) :: self !< Weights. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - endsubroutine compute_with_stencil_of_rank_1_interface + endsubroutine compute_stencil_rank_1_interface - pure subroutine compute_with_stencil_of_rank_2_interface(self, stencil) + pure subroutine compute_stencil_rank_2_interface(self, stencil, values) !< Compute beta. import :: weights_object, RPP - class(weights_object), intent(inout) :: self !< Weights. - real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. - endsubroutine compute_with_stencil_of_rank_2_interface + class(weights_object), intent(in) :: self !< Weights. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:,0:) !< Weights values of stencil interpolations. + endsubroutine compute_stencil_rank_2_interface pure subroutine smoothness_indicators_of_rank_1_interface(self, si) !< Return smoothness indicators. @@ -63,5 +63,4 @@ pure subroutine smoothness_indicators_of_rank_2_interface(self, si) real(RPP), intent(out) :: si(:,:) !< Smoothness indicators. endsubroutine smoothness_indicators_of_rank_2_interface endinterface - endmodule wenoof_weights_object diff --git a/src/lib/concrete_objects/wenoof_alpha_int_js.F90 b/src/lib/concrete_objects/wenoof_alpha_int_js.F90 index d59177d..be87620 100644 --- a/src/lib/concrete_objects/wenoof_alpha_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_alpha_int_js.F90 @@ -31,10 +31,11 @@ module wenoof_alpha_int_js !< 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_int !< Compute alpha. - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute_interpolate !< Compute alpha (interpolate). + procedure, pass(self) :: compute_reconstruct !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return alpha string-description. + procedure, pass(self) :: destroy !< Destroy alpha. endtype alpha_int_js contains @@ -53,21 +54,28 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_alpha_int(self, beta, kappa) + pure subroutine compute_interpolate(self, beta, kappa, values) !< Compute alpha. - class(alpha_int_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) :: s1 !< Counter. + class(alpha_int_js), intent(in) :: self !< Alpha coefficient. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + integer(I_P) :: s1 !< Counter. - associate(val => self%values_rank_1, val_sum => self%values_sum_rank_1) - val_sum = 0._RPP - do s1=0, self%S - 1 ! stencil loops - val(s1) = kappa%values_rank_1(s1)/(self%eps + beta%values_rank_1(s1)) ** self%S - val_sum = val_sum + val(s1) - enddo - endassociate - endsubroutine compute_alpha_int + do s1=0, self%S - 1 ! stencil loops + values(s1) = kappa(s1) / (self%eps + beta(s1)) ** self%S + enddo + endsubroutine compute_interpolate + + pure subroutine compute_reconstruct(self, beta, kappa, values) + !< Compute alpha. + class(alpha_int_js), intent(in) :: self !< Alpha coefficient. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + + ! Empty procedure. + endsubroutine compute_reconstruct pure function description(self) result(string) !< Return alpha string-descripition. diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 b/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 index 86f9ec4..a409691 100644 --- a/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 @@ -10,10 +10,10 @@ module wenoof_alpha_rec_js #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_alpha_object -use wenoof_base_object -use wenoof_beta_object -use wenoof_kappa_object +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_base_object, only : base_object_constructor +use wenoof_beta_object, only : beta_object +use wenoof_kappa_object, only : kappa_object implicit none private @@ -31,10 +31,11 @@ module wenoof_alpha_rec_js !< 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_rec !< Compute alpha. - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute_interpolate !< Compute alpha (interpolate). + procedure, pass(self) :: compute_reconstruct !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return alpha string-description. + procedure, pass(self) :: destroy !< Destroy alpha. endtype alpha_rec_js contains @@ -46,31 +47,32 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - allocate(self%values_sum_rank_2(1:2)) - associate(val => self%values_rank_2, val_sum => self%values_sum_rank_2) - val = 0._RPP - val_sum = 0._RPP - endassociate endsubroutine create - pure subroutine compute_alpha_rec(self, beta, kappa) + pure subroutine compute_interpolate_interface(self, beta, kappa, values) !< 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. + class(alpha_rec_js), intent(in) :: self !< Alpha coefficient. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. - associate(val => self%values_rank_2, val_sum => self%values_sum_rank_2) - val_sum = 0._RPP - do s1=0, self%S - 1 ! stencil loops - do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - val(f, s1) = kappa%values_rank_2(f, s1)/(self%eps + beta%values_rank_2(f, s1)) ** self%S - val_sum(f) = val_sum(f) + val(f, s1) - enddo + ! Empty procedure. + endsubroutine compute_interpolate_interface + + pure subroutine compute_reconstruct(self, beta, kappa, values) + !< Compute alpha. + class(alpha_rec_js), intent(in) :: self !< Alpha coefficient. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + integer(I_P) :: f, s1 !< Counters. + + do s1=0, self%S - 1 ! stencil loops + do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + values(f, s1) = kappa(f, s1) / (self%eps + beta(f, s1)) ** self%S enddo - endassociate - endsubroutine compute_alpha_rec + enddo + endsubroutine compute_reconstruct pure function description(self) result(string) !< Return alpha string-descripition. @@ -88,7 +90,5 @@ elemental subroutine destroy(self) class(alpha_rec_js), intent(inout) :: self !< Alpha. call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) - if (allocated(self%values_sum_rank_2)) deallocate(self%values_sum_rank_2) endsubroutine destroy endmodule wenoof_alpha_rec_js diff --git a/src/lib/concrete_objects/wenoof_beta_int_js.F90 b/src/lib/concrete_objects/wenoof_beta_int_js.F90 index 9c3ca18..1c1b215 100644 --- a/src/lib/concrete_objects/wenoof_beta_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_beta_int_js.F90 @@ -33,11 +33,11 @@ module wenoof_beta_int_js real(RPP), allocatable :: coef(:,:,:) !< Beta coefficients [0:S-1,0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create beta. - procedure, pass(self) :: compute_with_stencil_of_rank_1 !< Compute beta. - procedure, pass(self) :: compute_with_stencil_of_rank_2 !< Compute beta. - procedure, pass(self) :: description !< Return beta string-description. - procedure, pass(self) :: destroy !< Destroy beta. + procedure, pass(self) :: create !< Create beta. + procedure, pass(self) :: compute_stencil_rank_1 !< Compute beta. + procedure, pass(self) :: compute_stencil_rank_2 !< Compute beta. + procedure, pass(self) :: description !< Return beta string-description. + procedure, pass(self) :: destroy !< Destroy beta. endtype beta_int_js contains @@ -2372,7 +2372,7 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_with_stencil_of_rank_1(self, stencil) + pure subroutine compute_stencil_rank_1(self, stencil) !< Compute beta. class(beta_int_js), intent(inout) :: self !< Beta. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. @@ -2388,15 +2388,16 @@ pure subroutine compute_with_stencil_of_rank_1(self, stencil) enddo enddo endassociate - endsubroutine compute_with_stencil_of_rank_1 + endsubroutine compute_stencil_rank_1 - pure subroutine compute_with_stencil_of_rank_2(self, stencil) + pure subroutine compute_stencil_rank_2(self, stencil, values) !< Compute beta. - class(beta_int_js), intent(inout) :: self !< Beta. - real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + class(beta_int_js), intent(in) :: self !< Beta. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:,0:) !< Beta values [1:2,0:S-1]. ! Empty subroutine. - endsubroutine compute_with_stencil_of_rank_2 + endsubroutine compute_stencil_rank_2 pure function description(self) result(string) !< Return beta string-description. diff --git a/src/lib/concrete_objects/wenoof_beta_rec_js.F90 b/src/lib/concrete_objects/wenoof_beta_rec_js.F90 index 94a34e0..2a95068 100644 --- a/src/lib/concrete_objects/wenoof_beta_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_beta_rec_js.F90 @@ -35,24 +35,22 @@ module wenoof_beta_rec_js real(RPP), allocatable :: coef(:,:,:) !< Beta coefficients [0:S-1,0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create beta. - procedure, pass(self) :: compute_with_stencil_of_rank_1 !< Compute beta. - procedure, pass(self) :: compute_with_stencil_of_rank_2 !< Compute beta. - procedure, pass(self) :: description !< Return beta string-description. - procedure, pass(self) :: destroy !< Destroy beta. + procedure, pass(self) :: create !< Create beta. + procedure, pass(self) :: compute_stencil_rank_1 !< Compute beta. + procedure, pass(self) :: compute_stencil_rank_2 !< Compute beta. + procedure, pass(self) :: description !< Return beta string-description. + procedure, pass(self) :: destroy !< Destroy beta. endtype beta_rec_js contains ! 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. + class(beta_rec_js), intent(inout) :: self !< Beta. + class(base_object_constructor), intent(in) :: constructor !< Beta constructor. call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - self%values_rank_2 = 0._RPP allocate(self%coef(0:self%S - 1, 0:self%S - 1, 0:self%S - 1)) associate(c => self%coef) select case(self%S) @@ -2374,33 +2372,32 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_with_stencil_of_rank_1(self, stencil) + pure subroutine compute_stencil_rank_1(self, stencil) !< Compute beta. class(beta_rec_js), intent(inout) :: self !< Beta. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. ! Empty routine. - endsubroutine compute_with_stencil_of_rank_1 + endsubroutine compute_stencil_rank_1 - pure subroutine compute_with_stencil_of_rank_2(self, stencil) + pure subroutine compute_stencil_rank_2(self, stencil, values) !< Compute beta. - class(beta_rec_js), intent(inout) :: self !< Beta. - real(RPP), 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. - - associate(val => self%values_rank_2) - do s1=0, self%S - 1 ! stencils loop - do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - val(f, s1) = 0._RPP - do s2=0, self%S - 1 - do s3=0, self%S - 1 - val(f, s1) = val(f, s1) + self%coef(s3, s2, s1) * stencil(f, s1 - s3) * stencil(f, s1 - s2) - enddo + class(beta_rec_js), intent(in) :: self !< Beta. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:,0:) !< Beta values [1:2,0:S-1]. + integer(I_P) :: s1, s2, s3, f !< Counters. + + do s1=0, self%S - 1 ! stencils loop + do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + values(f, s1) = 0._RPP + do s2=0, self%S - 1 + do s3=0, self%S - 1 + values(f, s1) = val(f, s1) + self%coef(s3, s2, s1) * stencil(f, s1 - s3) * stencil(f, s1 - s2) enddo enddo enddo - endassociate - endsubroutine compute_with_stencil_of_rank_2 + enddo + endsubroutine compute_stencil_rank_2 pure function description(self) result(string) !< Return beta string-description. @@ -2418,7 +2415,6 @@ elemental subroutine destroy(self) class(beta_rec_js), intent(inout) :: self !< Beta. call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) if (allocated(self%coef)) deallocate(self%coef) endsubroutine destroy endmodule wenoof_beta_rec_js diff --git a/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 b/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 index ae072e2..a1b7c02 100644 --- a/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 @@ -32,11 +32,11 @@ module wenoof_interpolations_int_js real(RPP), allocatable :: coef(:,:) !< Polynomial coefficients [0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create interpolations. - procedure, pass(self) :: compute_with_stencil_of_rank_1 !< Compute interpolations. - procedure, pass(self) :: compute_with_stencil_of_rank_2 !< Compute interpolations. - procedure, pass(self) :: description !< Return interpolations string-description. - procedure, pass(self) :: destroy !< Destroy interpolations. + procedure, pass(self) :: create !< Create interpolations. + procedure, pass(self) :: compute_stencil_rank_1 !< Compute interpolations. + procedure, pass(self) :: compute_stencil_rank_2 !< Compute interpolations. + procedure, pass(self) :: description !< Return interpolations string-description. + procedure, pass(self) :: destroy !< Destroy interpolations. endtype interpolations_int_js contains @@ -346,7 +346,7 @@ subroutine create(self, constructor) endselect endsubroutine create - pure subroutine compute_with_stencil_of_rank_1(self, stencil) + pure subroutine compute_stencil_rank_1(self, stencil) !< Compute interpolations. class(interpolations_int_js), intent(inout) :: self !< Interpolations. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. @@ -361,15 +361,16 @@ pure subroutine compute_with_stencil_of_rank_1(self, stencil) enddo enddo endassociate - endsubroutine compute_with_stencil_of_rank_1 + endsubroutine compute_stencil_rank_1 - pure subroutine compute_with_stencil_of_rank_2(self, stencil) + pure subroutine compute_stencil_rank_2(self, stencil, values) !< Compute interpolations. - class(interpolations_int_js), intent(inout) :: self !< Interpolations. - real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + class(interpolations_int_js), intent(in) :: self !< Interpolations. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:, 0:) !< Interpolations values. ! Empty Subroutine. - endsubroutine compute_with_stencil_of_rank_2 + endsubroutine compute_stencil_rank_2 pure function description(self) result(string) !< Return interpolations string-description. diff --git a/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 b/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 index 105792d..b4fcd44 100644 --- a/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 @@ -35,11 +35,11 @@ module wenoof_interpolations_rec_js real(RPP), allocatable :: coef(:,:,:) !< Polynomial coefficients [1:2,0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create interpolations. - procedure, pass(self) :: compute_with_stencil_of_rank_1 !< Compute interpolations. - procedure, pass(self) :: compute_with_stencil_of_rank_2 !< Compute interpolations. - procedure, pass(self) :: description !< Return interpolations string-description. - procedure, pass(self) :: destroy !< Destroy interpolations. + procedure, pass(self) :: create !< Create interpolations. + procedure, pass(self) :: compute_stencil_rank_1 !< Compute interpolations. + procedure, pass(self) :: compute_stencil_rank_2 !< Compute interpolations. + procedure, pass(self) :: description !< Return interpolations string-description. + procedure, pass(self) :: destroy !< Destroy interpolations. endtype interpolations_rec_js contains @@ -51,8 +51,6 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - self%values_rank_2 = 0._RPP allocate(self%coef(1:2, 0:self%S - 1, 0:self%S - 1)) associate(c => self%coef) select case(self%S) @@ -451,33 +449,32 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_with_stencil_of_rank_1(self, stencil) + pure subroutine compute_stencil_rank_1(self, stencil) !< Compute interpolations. class(interpolations_rec_js), intent(inout) :: self !< Interpolations. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. ! Empty Subroutine. - endsubroutine compute_with_stencil_of_rank_1 + endsubroutine compute_stencil_rank_1 - pure subroutine compute_with_stencil_of_rank_2(self, stencil) + pure subroutine compute_stencil_rank_2(self, stencil, values) !< Compute interpolations. - class(interpolations_rec_js), intent(inout) :: self !< Interpolations. - real(RPP), 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. + class(interpolations_rec_js), intent(in) :: self !< Interpolations. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:, 0:) !< Interpolations values. + integer(I_P) :: s1 !< Counter. + integer(I_P) :: s2 !< Counter. + integer(I_P) :: f !< Counter. - associate(val => self%values_rank_2) - val = 0._RPP + values = 0._RPP do s1=0, self%S - 1 ! stencils loop do s2=0, self%S - 1 ! values loop do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - val(f, s1) = val(f, s1) + self%coef(f, s2, s1) * stencil(f, -s2 + s1) + values(f, s1) = values(f, s1) + self%coef(f, s2, s1) * stencil(f, -s2 + s1) enddo enddo enddo - endassociate - endsubroutine compute_with_stencil_of_rank_2 + endsubroutine compute_stencil_rank_2 pure function description(self) result(string) !< Return interpolations string-description. @@ -495,7 +492,6 @@ elemental subroutine destroy(self) class(interpolations_rec_js), intent(inout) :: self !< Interpolations. call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) if (allocated(self%coef)) deallocate(self%coef) endsubroutine destroy 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 index 05d3b1e..7aac609 100644 --- a/src/lib/concrete_objects/wenoof_interpolator_js.F90 +++ b/src/lib/concrete_objects/wenoof_interpolator_js.F90 @@ -34,13 +34,13 @@ module wenoof_interpolator_js !< 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 interpolator. - procedure, pass(self) :: description !< Return interpolator string-description. - procedure, pass(self) :: destroy !< Destroy interpolator. - procedure, pass(self) :: interpolate_with_stencil_of_rank_1_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_with_stencil_of_rank_2_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_with_stencil_of_rank_1_debug !< Interpolate values (providing also debug values). - procedure, pass(self) :: interpolate_with_stencil_of_rank_2_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: create !< Create interpolator. + procedure, pass(self) :: description !< Return interpolator string-description. + procedure, pass(self) :: destroy !< Destroy interpolator. + procedure, pass(self) :: interpolate_stencil_rank_1_standard !< Interpolate values (without providing debug values). + procedure, pass(self) :: interpolate_stencil_rank_2_standard !< Interpolate values (without providing debug values). + procedure, pass(self) :: interpolate_stencil_rank_1_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: interpolate_stencil_rank_2_debug !< Interpolate values (providing also debug values). endtype interpolator_js contains @@ -81,7 +81,7 @@ elemental subroutine destroy(self) if (allocated(self%weights)) deallocate(self%weights) endsubroutine destroy - pure subroutine interpolate_with_stencil_of_rank_1_debug(self, stencil, interpolation, si, weights) + pure subroutine interpolate_stencil_rank_1_debug(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). class(interpolator_js), intent(inout) :: self !< Interpolator. real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. @@ -92,9 +92,9 @@ pure subroutine interpolate_with_stencil_of_rank_1_debug(self, stencil, interpol call self%interpolate(stencil=stencil, interpolation=interpolation) call self%weights%smoothness_indicators_of_rank_1(si=si) weights = self%weights%values_rank_1 - endsubroutine interpolate_with_stencil_of_rank_1_debug + endsubroutine interpolate_stencil_rank_1_debug - pure subroutine interpolate_with_stencil_of_rank_2_debug(self, stencil, interpolation, si, weights) + pure subroutine interpolate_stencil_rank_2_debug(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). class(interpolator_js), intent(inout) :: self !< Reconstructor. real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. @@ -103,9 +103,9 @@ pure subroutine interpolate_with_stencil_of_rank_2_debug(self, stencil, interpol real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. ! Empty subroutine. - endsubroutine interpolate_with_stencil_of_rank_2_debug + endsubroutine interpolate_stencil_rank_2_debug - pure subroutine interpolate_with_stencil_of_rank_1_standard(self, stencil, interpolation) + pure subroutine interpolate_stencil_rank_1_standard(self, stencil, interpolation) !< Interpolate values (without providing debug values). class(interpolator_js), intent(inout) :: self !< Interpolator. real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. @@ -118,14 +118,14 @@ pure subroutine interpolate_with_stencil_of_rank_1_standard(self, stencil, inter do s=0, self%S - 1 ! stencils loop interpolation = interpolation + self%weights%values_rank_1(s) * self%interpolations%values_rank_1(s) enddo - endsubroutine interpolate_with_stencil_of_rank_1_standard + endsubroutine interpolate_stencil_rank_1_standard - pure subroutine interpolate_with_stencil_of_rank_2_standard(self, stencil, interpolation) + pure subroutine interpolate_stencil_rank_2_standard(self, stencil, interpolation) !< Interpolate values (without providing debug values). class(interpolator_js), intent(inout) :: self !< Reconstructor. real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. ! Empty subroutine. - endsubroutine interpolate_with_stencil_of_rank_2_standard + endsubroutine interpolate_stencil_rank_2_standard endmodule wenoof_interpolator_js diff --git a/src/lib/concrete_objects/wenoof_kappa_int_js.F90 b/src/lib/concrete_objects/wenoof_kappa_int_js.F90 index 8226436..4b5705c 100644 --- a/src/lib/concrete_objects/wenoof_kappa_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_kappa_int_js.F90 @@ -70,7 +70,7 @@ pure subroutine compute_kappa_rec(self) !< Compute kappa. class(kappa_int_js), intent(inout) :: self !< Kappa. - ! Empty subroutine + ! Empty procedure. endsubroutine compute_kappa_rec pure subroutine compute_kappa_int(self, stencil, x_target) diff --git a/src/lib/concrete_objects/wenoof_reconstructor_js.F90 b/src/lib/concrete_objects/wenoof_reconstructor_js.F90 index 171cfdd..af182b6 100644 --- a/src/lib/concrete_objects/wenoof_reconstructor_js.F90 +++ b/src/lib/concrete_objects/wenoof_reconstructor_js.F90 @@ -34,13 +34,13 @@ module wenoof_reconstructor_js !< 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_with_stencil_of_rank_1_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_with_stencil_of_rank_2_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_with_stencil_of_rank_1_debug !< Interpolate values (providing also debug values). - procedure, pass(self) :: interpolate_with_stencil_of_rank_2_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: create !< Create reconstructor. + procedure, pass(self) :: description !< Return reconstructor string-description. + procedure, pass(self) :: destroy !< Destroy reconstructor. + procedure, pass(self) :: interpolate_stencil_rank_1_standard !< Interpolate values (without providing debug values). + procedure, pass(self) :: interpolate_stencil_rank_2_standard !< Interpolate values (without providing debug values). + procedure, pass(self) :: interpolate_stencil_rank_1_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: interpolate_stencil_rank_2_debug !< Interpolate values (providing also debug values). endtype reconstructor_js contains @@ -81,7 +81,7 @@ elemental subroutine destroy(self) if (allocated(self%weights)) deallocate(self%weights) endsubroutine destroy - pure subroutine interpolate_with_stencil_of_rank_1_debug(self, stencil, interpolation, si, weights) + pure subroutine interpolate_stencil_rank_1_debug(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). class(reconstructor_js), intent(inout) :: self !< Reconstructor. real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. @@ -90,9 +90,9 @@ pure subroutine interpolate_with_stencil_of_rank_1_debug(self, stencil, interpol real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [0:S-1]. ! Empty subroutine. - endsubroutine interpolate_with_stencil_of_rank_1_debug + endsubroutine interpolate_stencil_rank_1_debug - pure subroutine interpolate_with_stencil_of_rank_2_debug(self, stencil, interpolation, si, weights) + pure subroutine interpolate_stencil_rank_2_debug(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). class(reconstructor_js), intent(inout) :: self !< Reconstructor. real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. @@ -103,32 +103,35 @@ pure subroutine interpolate_with_stencil_of_rank_2_debug(self, stencil, interpol call self%interpolate(stencil=stencil, interpolation=interpolation) call self%weights%smoothness_indicators_of_rank_2(si=si) !si = self%weights%smoothness_indicators() - weights = self%weights%values_rank_2 - endsubroutine interpolate_with_stencil_of_rank_2_debug + ! weights = self%weights%values_rank_2 + endsubroutine interpolate_stencil_rank_2_debug - pure subroutine interpolate_with_stencil_of_rank_1_standard(self, stencil, interpolation) + pure subroutine interpolate_stencil_rank_1_standard(self, stencil, interpolation) !< Interpolate values (without providing debug values). class(reconstructor_js), intent(inout) :: self !< Reconstructor. real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. real(RPP), intent(out) :: interpolation !< Result of the interpolation. ! Empty subroutine. - endsubroutine interpolate_with_stencil_of_rank_1_standard + endsubroutine interpolate_stencil_rank_1_standard - pure subroutine interpolate_with_stencil_of_rank_2_standard(self, stencil, interpolation) + pure subroutine interpolate_stencil_rank_2_standard(self, stencil, interpolation) !< Interpolate values (without providing debug values). - class(reconstructor_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), 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) + class(reconstructor_js), intent(inout) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(RPP) :: interpolations(1:2, 0:self%S - 1) !< Stencils interpolations. + real(RPP) :: weights(1:2, 0:self%S - 1) !< Weights of stencils interpolations. + integer(I_P) :: f, s !< Counters. + + call self%interpolations%compute(stencil=stencil, values=interpolations) + call self%weights%compute(stencil=stencil, values=weights) interpolation = 0._RPP do s=0, self%S - 1 ! stencils loop do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - interpolation(f) = interpolation(f) + self%weights%values_rank_2(f, s) * self%interpolations%values_rank_2(f, s) + ! interpolation(f) = interpolation(f) + self%weights%values_rank_2(f, s) * self%interpolations%values_rank_2(f, s) + interpolation(f) = interpolation(f) + weights(f, s) * interpolations(f, s) enddo enddo - endsubroutine interpolate_with_stencil_of_rank_2_standard + endsubroutine interpolate_stencil_rank_2_standard endmodule wenoof_reconstructor_js diff --git a/src/lib/concrete_objects/wenoof_weights_int_js.F90 b/src/lib/concrete_objects/wenoof_weights_int_js.F90 index cc4a7d2..2de2047 100644 --- a/src/lib/concrete_objects/wenoof_weights_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_weights_int_js.F90 @@ -45,14 +45,14 @@ module wenoof_weights_int_js !< 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). + 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_with_stencil_of_rank_1 !< Compute weights. - procedure, pass(self) :: compute_with_stencil_of_rank_2 !< Compute weights. + procedure, pass(self) :: compute_stencil_rank_1 !< Compute weights. + procedure, pass(self) :: compute_stencil_rank_2 !< Compute weights. procedure, pass(self) :: description !< Return weights string-description. procedure, pass(self) :: destroy !< Destroy weights. procedure, pass(self) :: smoothness_indicators_of_rank_1 !< Return smoothness indicators. @@ -106,7 +106,7 @@ subroutine create(self, constructor) endselect endsubroutine create - pure subroutine compute_with_stencil_of_rank_1(self, stencil) + pure subroutine compute_stencil_rank_1(self, stencil) !< Compute weights. class(weights_int_js), intent(inout) :: self !< Weights. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. @@ -117,15 +117,16 @@ pure subroutine compute_with_stencil_of_rank_1(self, stencil) do s=0, self%S - 1 ! stencils loop self%values_rank_1(s) = self%alpha%values_rank_1(s) / self%alpha%values_sum_rank_1 enddo - endsubroutine compute_with_stencil_of_rank_1 + endsubroutine compute_stencil_rank_1 - pure subroutine compute_with_stencil_of_rank_2(self, stencil) + pure subroutine compute_stencil_rank_2(self, stencil, values) !< Compute weights. - class(weights_int_js), intent(inout) :: self !< Weights. - real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + class(weights_int_js), intent(in) :: self !< Weights. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:,0:) !< Weights values of stencil interpolations. ! Empty routine. - endsubroutine compute_with_stencil_of_rank_2 + endsubroutine compute_stencil_rank_2 pure function description(self) result(string) !< Return string-description of weights. diff --git a/src/lib/concrete_objects/wenoof_weights_rec_js.F90 b/src/lib/concrete_objects/wenoof_weights_rec_js.F90 index 5f5056e..7191b50 100644 --- a/src/lib/concrete_objects/wenoof_weights_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_weights_rec_js.F90 @@ -45,14 +45,14 @@ module wenoof_weights_rec_js !< 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). + 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_with_stencil_of_rank_1 !< Compute weights. - procedure, pass(self) :: compute_with_stencil_of_rank_2 !< Compute weights. + procedure, pass(self) :: compute_stencil_rank_1 !< Compute weights. + procedure, pass(self) :: compute_stencil_rank_2 !< Compute weights. procedure, pass(self) :: description !< Return weights string-description. procedure, pass(self) :: destroy !< Destroy weights. procedure, pass(self) :: smoothness_indicators_of_rank_1 !< Return smoothness indicators. @@ -71,8 +71,6 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - self%values_rank_2 = 0._RPP select type(constructor) type is(weights_rec_js_constructor) associate(alpha_constructor=>constructor%alpha_constructor, & @@ -106,28 +104,30 @@ subroutine create(self, constructor) endselect endsubroutine create - pure subroutine compute_with_stencil_of_rank_1(self, stencil) + pure subroutine compute_stencil_rank_1(self, stencil) !< Compute weights. class(weights_rec_js), intent(inout) :: self !< Weights. real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. ! Empty routine. - endsubroutine compute_with_stencil_of_rank_1 + endsubroutine compute_stencil_rank_1 - pure subroutine compute_with_stencil_of_rank_2(self, stencil) + pure subroutine compute_stencil_rank_2(self, stencil, values) !< Compute weights. - class(weights_rec_js), intent(inout) :: self !< Weights. - real(RPP), 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) + class(weights_rec_js), intent(in) :: self !< Weights. + real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: values(1:,0:) !< Weights values of stencil interpolations. + real(RPP) :: beta(1:2,0:self%S-1) !< Beta values. + integer(I_P) :: f, s !< Counters. + + call self%beta%compute(stencil=stencil, values=beta) + ! call self%alpha%compute(beta=self%beta, kappa=self%kappa) do s=0, self%S - 1 ! stencils loop do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - self%values_rank_2(f, s) = self%alpha%values_rank_2(f, s) / self%alpha%values_sum_rank_2(f) + values(f, s) = self%alpha%values_rank_2(f, s) / self%alpha%values_sum_rank_2(f) enddo enddo - endsubroutine compute_with_stencil_of_rank_2 + endsubroutine compute_stencil_rank_2 pure function description(self) result(string) !< Return string-description of weights. @@ -145,7 +145,6 @@ elemental subroutine destroy(self) class(weights_rec_js), intent(inout) :: self !< Weights. call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) if (allocated(self%alpha)) deallocate(self%alpha) if (allocated(self%beta)) deallocate(self%beta) if (allocated(self%kappa)) deallocate(self%kappa) diff --git a/src/lib/wenoof.F90 b/src/lib/wenoof.F90 index 6484a4c..c806d83 100644 --- a/src/lib/wenoof.F90 +++ b/src/lib/wenoof.F90 @@ -7,8 +7,8 @@ module wenoof #else use penf, only: I_P, RPP=>R8P #endif -use wenoof_interpolator_object -use wenoof_objects_factory +use wenoof_interpolator_object, only : interpolator_object +use wenoof_objects_factory, only : objects_factory implicit none private @@ -50,7 +50,7 @@ subroutine wenoof_create_interpolator(interpolator_type, S, x_target, interpolat error stop 'error: x_target must be between -0.5 and 0.5, that represent left and right cell interfaces' endif allocate(stencil(1-S:S-1)) - do i=0,2*S-2 + do i=0, 2 * S - 2 stencil(-S+1+i) = 1.0_RPP - S + i enddo call factory%create(interpolator_type=interpolator_type, & From e7a79c20b9f27bdf559584d0e49f59c8365946de Mon Sep 17 00:00:00 2001 From: Stefano Zaghi Date: Mon, 10 Apr 2017 17:25:10 +0200 Subject: [PATCH 3/3] API refactoring API refactoring The API has been refactored: all objects (alpha, beta, kappa, interpolations and weights) now compute their values into volatile buffers thus the interpolator can be defined as `intent(in)` when used for computing an interpolation. Why: Efficiency reasons: now an interpolator can be defined as a member of user derived type and used to interpolate values even into procedure where the user derived type is defined as `intent(in)` without the unnecessary re-initialization or unsafe *global* definition. This change addresses the need by: API of all objects has been changed. Side effects: Backward compatibility lost. --- quarantena/wenoof_alpha_int_m.F90 | 131 -------- quarantena/wenoof_alpha_rec_m.F90 | 136 -------- .../abstract_objects/wenoof_alpha_object.F90 | 31 +- .../abstract_objects/wenoof_base_object.F90 | 7 +- .../abstract_objects/wenoof_beta_object.F90 | 26 +- .../wenoof_interpolations_object.F90 | 31 +- .../wenoof_interpolator_object.F90 | 75 +++-- .../abstract_objects/wenoof_kappa_object.F90 | 39 ++- .../wenoof_weights_object.F90 | 48 +-- .../concrete_objects/wenoof_alpha_int_js.F90 | 56 ++-- .../concrete_objects/wenoof_alpha_int_m.F90 | 132 ++++++++ .../concrete_objects}/wenoof_alpha_int_z.F90 | 85 +++-- .../concrete_objects/wenoof_alpha_rec_js.F90 | 46 +-- .../concrete_objects/wenoof_alpha_rec_m.F90 | 135 ++++++++ .../concrete_objects}/wenoof_alpha_rec_z.F90 | 87 +++-- .../concrete_objects/wenoof_beta_int_js.F90 | 81 +++-- .../concrete_objects/wenoof_beta_rec_js.F90 | 62 ++-- .../wenoof_interpolations_int_js.F90 | 71 ++--- .../wenoof_interpolations_rec_js.F90 | 56 ++-- .../wenoof_interpolator_js.F90 | 141 +++++---- .../concrete_objects/wenoof_kappa_int_js.F90 | 299 +++++++++--------- .../concrete_objects/wenoof_kappa_rec_js.F90 | 298 ++++++++--------- .../wenoof_reconstructor_js.F90 | 132 ++++---- .../wenoof_weights_int_js.F90 | 138 ++++---- .../wenoof_weights_rec_js.F90 | 133 ++++---- ...a_factory.f90 => wenoof_alpha_factory.F90} | 16 +- src/lib/factories/wenoof_beta_factory.f90 | 12 +- ....f90 => wenoof_interpolations_factory.F90} | 6 +- .../factories/wenoof_interpolator_factory.f90 | 40 +-- src/lib/factories/wenoof_kappa_factory.f90 | 23 +- src/lib/factories/wenoof_objects_factory.f90 | 24 +- src/lib/factories/wenoof_weights_factory.f90 | 16 +- 32 files changed, 1280 insertions(+), 1333 deletions(-) delete mode 100644 quarantena/wenoof_alpha_int_m.F90 delete mode 100644 quarantena/wenoof_alpha_rec_m.F90 create mode 100644 src/lib/concrete_objects/wenoof_alpha_int_m.F90 rename {quarantena => src/lib/concrete_objects}/wenoof_alpha_int_z.F90 (56%) create mode 100644 src/lib/concrete_objects/wenoof_alpha_rec_m.F90 rename {quarantena => src/lib/concrete_objects}/wenoof_alpha_rec_z.F90 (56%) rename src/lib/factories/{wenoof_alpha_factory.f90 => wenoof_alpha_factory.F90} (87%) rename src/lib/factories/{wenoof_interpolations_factory.f90 => wenoof_interpolations_factory.F90} (91%) diff --git a/quarantena/wenoof_alpha_int_m.F90 b/quarantena/wenoof_alpha_int_m.F90 deleted file mode 100644 index 40eb910..0000000 --- a/quarantena/wenoof_alpha_int_m.F90 +++ /dev/null @@ -1,131 +0,0 @@ -!< Henrick alpha (non linear weights) object. -module wenoof_alpha_int_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 - -#ifdef r16p -use penf, only: I_P, RPP=>R16P, str -#else -use penf, only: I_P, RPP=>R8P, str -#endif -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_int_m -public :: alpha_int_m_constructor - -type, extends(alpha_object_constructor) :: alpha_int_m_constructor - !< Henrick alpha (non linear weights) object constructor. - character(len=:), allocatable :: base_type !< Base alpha coefficient type. -endtype alpha_int_m_constructor - -type, extends(alpha_object) :: alpha_int_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. - real(RPP), allocatable :: values(:) !< Alpha coefficients [0:S-1]. - real(RPP) :: values_sum !< Sum of alpha coefficients. - 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_int !< Compute alpha. - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. -endtype alpha_int_m - -contains - ! deferred public methods - subroutine create(self, constructor) - !< Create alpha. - class(alpha_int_m), intent(inout) :: self !< Alpha. - class(base_object_constructor), intent(in) :: constructor !< Alpha constructor. - - call self%destroy - call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - associate(val => self%values_rank_1, val_sum => self%values_sum_rank_1) - val = 0._RPP - val_sum = 0._RPP - endassociate - select type(constructor) - type is(alpha_int_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_alpha_int(self, beta, kappa) - !< Compute alpha. - class(alpha_int_m), intent(inout) :: self !< Alpha. - class(beta_object), intent(in) :: beta !< Beta. - class(kappa_object), intent(in) :: kappa !< Kappa. - real(RPP) :: kappa_base !< Kappa evaluated from the base alphas. - integer(I_P) :: s1 !< Counter. - - associate(val => self%values_rank_1, val_sum => self%values_sum_rank_1) - val_sum = 0._RPP - call self%alpha_base%compute(beta=beta, kappa=kappa) - do s1=0, self%S - 1 ! stencil loops - kappa_base = self%alpha_base%values_rank_1(s1) / self%alpha_base%values_sum_rank_1 - val(s1) = & - (kappa_base * (kappa%values_rank_1(s1) + kappa%values_rank_1(s1) * kappa%values_rank_1(s1) - & - 3._RPP * kappa%values_rank_1(s1) * kappa_base + kappa_base * kappa_base)) / & - (kappa%values_rank_1(s1) * kappa%values_rank_1(s1) + kappa_base * & - (1._RPP - 2._RPP * kappa%values_rank_1(s1))) - val_sum = val_sum + val(s1) - enddo - endassociate - endsubroutine compute_alpha_int - - pure function description(self) result(string) - !< Return alpha string-descripition. - class(alpha_int_m), intent(in) :: self !< Alpha. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Henrick alpha coefficients for reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//' - eps = '//trim(str(self%eps))//nl - associate(alpha_base=>self%alpha_base) - select type(alpha_base) - type is(alpha_rec_js) - string = string//' - base-mapped-alpha type = Jiang-Shu' - type is(alpha_rec_z) - string = string//' - base-mapped-alpha type = Bogeg' - endselect - endassociate - endfunction description - - elemental subroutine destroy(self) - !< Destroy alpha. - class(alpha_int_m), intent(inout) :: self !< Alpha. - - call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) - if (allocated(self%alpha_base)) deallocate(self%alpha_base) - endsubroutine destroy -endmodule wenoof_alpha_int_m diff --git a/quarantena/wenoof_alpha_rec_m.F90 b/quarantena/wenoof_alpha_rec_m.F90 deleted file mode 100644 index 7071b3e..0000000 --- a/quarantena/wenoof_alpha_rec_m.F90 +++ /dev/null @@ -1,136 +0,0 @@ -!< 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 - -#ifdef r16p -use penf, only: I_P, RPP=>R16P, str -#else -use penf, only: I_P, RPP=>R8P, str -#endif -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. - real(RPP), allocatable :: values(:,:) !< Alpha coefficients [1:2,0:S-1]. - real(RPP), allocatable :: values_sum(:) !< Sum of alpha coefficients [1:2]. - 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_rec !< 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) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - allocate(self%values_sum_rank_2(1:2)) - associate(val => self%values_rank_2, val_sum => self%values_sum_rank_2) - val = 0._RPP - val_sum = 0._RPP - endassociate - 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_alpha_rec(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. - real(RPP) :: kappa_base !< Kappa evaluated from the base alphas. - integer(I_P) :: f, s1 !< Counters. - - associate(val => self%values_rank_2, val_sum => self%values_sum_rank_2) - val_sum = 0._RPP - call self%alpha_base%compute(beta=beta, kappa=kappa) - do s1=0, self%S - 1 ! stencil loops - do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - kappa_base = self%alpha_base%values_rank_2(f, s1) / self%alpha_base%values_sum_rank_2(f) - val(f, s1) = & - (kappa_base * (kappa%values_rank_2(f, s1) + kappa%values_rank_2(f, s1) * kappa%values_rank_2(f, s1) - & - 3._RPP * kappa%values_rank_2(f, s1) * kappa_base + kappa_base * & - kappa_base)) / & - (kappa%values_rank_2(f, s1) * kappa%values_rank_2(f, s1) + kappa_base * & - (1._RPP - 2._RPP * kappa%values_rank_2(f, s1))) - val_sum(f) = val_sum(f) + val(f, s1) - enddo - enddo - endassociate - endsubroutine compute_alpha_rec - - pure function description(self) result(string) - !< Return alpha string-descripition. - class(alpha_rec_m), intent(in) :: self !< Alpha. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Henrick alpha coefficients for reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//' - eps = '//trim(str(self%eps))//nl - associate(alpha_base=>self%alpha_base) - select type(alpha_base) - type is(alpha_rec_js) - string = string//' - base-mapped-alpha type = Jiang-Shu' - type is(alpha_rec_z) - string = string//' - base-mapped-alpha type = Bogeg' - endselect - endassociate - endfunction description - - elemental subroutine destroy(self) - !< Destroy alpha. - class(alpha_rec_m), intent(inout) :: self !< Alpha. - - call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) - if (allocated(self%values_sum_rank_2)) deallocate(self%values_sum_rank_2) - if (allocated(self%alpha_base)) deallocate(self%alpha_base) - endsubroutine destroy -endmodule wenoof_alpha_rec_m diff --git a/src/lib/abstract_objects/wenoof_alpha_object.F90 b/src/lib/abstract_objects/wenoof_alpha_object.F90 index da24594..877b6d0 100644 --- a/src/lib/abstract_objects/wenoof_alpha_object.F90 +++ b/src/lib/abstract_objects/wenoof_alpha_object.F90 @@ -7,9 +7,7 @@ module wenoof_alpha_object #else use penf, only: RPP=>R8P #endif -use wenoof_base_object -use wenoof_beta_object -use wenoof_kappa_object +use wenoof_base_object, only : base_object, base_object_constructor implicit none private @@ -23,35 +21,32 @@ module wenoof_alpha_object type, extends(base_object), abstract :: alpha_object !< Abstract alpha (non linear weights) object. - ! real(RPP), allocatable :: values_rank_1(:) !< Alpha values [0:S-1]. - ! real(RPP) :: values_sum_rank_1 !< Sum of alpha coefficients. - ! real(RPP), allocatable :: values_rank_2(:,:) !< Alpha values [1:2,0:S-1]. - ! real(RPP), allocatable :: values_sum_rank_2(:) !< Sum of alpha coefficients [1:2]. contains + ! public methods + generic :: compute => compute_int, compute_rec !< Compute alpha. ! public deferred methods - procedure(compute_interpolate_interface), pass(self), deferred :: compute_interpolate !< Compute alpha (interpolate). - procedure(compute_reconstruct_interface), pass(self), deferred :: compute_reconstruct !< Compute alpha (reconstruct). - generic :: compute => compute_interpolate, compute_reconstruct + procedure(compute_int_interface), pass(self), deferred :: compute_int !< Compute alpha (interpolate). + procedure(compute_rec_interface), pass(self), deferred :: compute_rec !< Compute alpha (reconstruct). endtype alpha_object abstract interface !< Abstract interfaces of [[alpha_object]]. - pure subroutine compute_interpolate_interface(self, beta, kappa, values) - !< Compute alpha. - import :: alpha_object, beta_object, kappa_object + pure subroutine compute_int_interface(self, beta, kappa, values) + !< Compute alpha (interpolate). + import :: alpha_object, RPP class(alpha_object), intent(in) :: self !< Alpha. real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. - endsubroutine compute_interpolate_interface + endsubroutine compute_int_interface - pure subroutine compute_reconstruct_interface(self, beta, kappa, values) - !< Compute alpha. - import :: alpha_object, beta_object, kappa_object + pure subroutine compute_rec_interface(self, beta, kappa, values) + !< Compute alpha (reconstruct). + import :: alpha_object, RPP class(alpha_object), intent(in) :: self !< Alpha. real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. - endsubroutine compute_reconstruct_interface + endsubroutine compute_rec_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 index 50cbf9a..c0f31fa 100644 --- a/src/lib/abstract_objects/wenoof_base_object.F90 +++ b/src/lib/abstract_objects/wenoof_base_object.F90 @@ -52,11 +52,12 @@ subroutine create_interface(self, constructor) class(base_object_constructor), intent(in) :: constructor !< Object constructor. endsubroutine create_interface - pure function description_interface(self) result(string) + pure function description_interface(self, prefix) result(string) !< Return object string-description. import :: base_object - class(base_object), intent(in) :: self !< Object. - character(len=:), allocatable :: string !< String-description. + class(base_object), intent(in) :: self !< Object. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. endfunction description_interface elemental subroutine destroy_interface(self) diff --git a/src/lib/abstract_objects/wenoof_beta_object.F90 b/src/lib/abstract_objects/wenoof_beta_object.F90 index f15a9ad..a8e563d 100644 --- a/src/lib/abstract_objects/wenoof_beta_object.F90 +++ b/src/lib/abstract_objects/wenoof_beta_object.F90 @@ -7,7 +7,7 @@ module wenoof_beta_object #else use penf, only: RPP=>R8P #endif -use wenoof_base_object +use wenoof_base_object, only : base_object, base_object_constructor implicit none private @@ -20,30 +20,30 @@ module wenoof_beta_object type, extends(base_object), abstract :: beta_object !< Abstract Beta coefficients (smoothness indicators of stencil interpolations) object. - real(RPP), allocatable :: values_rank_1(:) !< Beta values [0:S-1]. contains ! public methods - generic :: compute => compute_stencil_rank_1, compute_stencil_rank_2 + generic :: compute => compute_int, compute_rec !< Compute beta. ! deferred public methods - procedure(compute_stencil_rank_1_interface), pass(self), deferred :: compute_stencil_rank_1 !< Compute beta. - procedure(compute_stencil_rank_2_interface), pass(self), deferred :: compute_stencil_rank_2 !< Compute beta. + procedure(compute_int_interface), pass(self), deferred :: compute_int !< Compute beta (interpolate). + procedure(compute_rec_interface), pass(self), deferred :: compute_rec !< Compute beta (reconstruct). endtype beta_object abstract interface !< Abstract interfaces of [[beta_object]]. - pure subroutine compute_stencil_rank_1_interface(self, stencil) - !< Compute beta. + pure subroutine compute_int_interface(self, stencil, values) + !< Compute beta (interpolate). import :: beta_object, RPP - class(beta_object), intent(inout) :: self !< Beta. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - endsubroutine compute_stencil_rank_1_interface + class(beta_object), intent(in) :: self !< Beta. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Beta values [0:S-1]. + endsubroutine compute_int_interface - pure subroutine compute_stencil_rank_2_interface(self, stencil, values) - !< Compute beta. + pure subroutine compute_rec_interface(self, stencil, values) + !< Compute beta (reconstruct). import :: beta_object, RPP class(beta_object), intent(in) :: self !< Beta. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:,0:) !< Beta values [1:2,0:S-1]. - endsubroutine compute_stencil_rank_2_interface + endsubroutine compute_rec_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 index cbe349f..4da73a2 100644 --- a/src/lib/abstract_objects/wenoof_interpolations_object.F90 +++ b/src/lib/abstract_objects/wenoof_interpolations_object.F90 @@ -7,7 +7,7 @@ module wenoof_interpolations_object #else use penf, only: RPP=>R8P #endif -use wenoof_base_object +use wenoof_base_object, only : base_object, base_object_constructor implicit none private @@ -16,37 +16,36 @@ module wenoof_interpolations_object type, extends(base_object_constructor), abstract :: interpolations_object_constructor !< Abstract interpolations object constructor. - real(RPP), allocatable :: stencil(:) !< Stencil used for interpolation, [1-S:S-1]. - real(RPP) :: x_target !< Coordinate of the interpolation point. + real(RPP), allocatable :: stencil(:) !< Stencil used for interpolation, [1-S:S-1]. + real(RPP) :: x_target !< Coordinate of the interpolation point. endtype interpolations_object_constructor type, extends(base_object), abstract :: interpolations_object !< Abstract interpolations object. - real(RPP), allocatable :: values_rank_1(:) !< Stencil interpolations values [0:S-1]. - ! real(RPP), allocatable :: values_rank_2(:,:) !< Stencil interpolations values [1:2,0:S-1]. contains ! public methods - generic :: compute => compute_stencil_rank_1, compute_stencil_rank_2 + generic :: compute => compute_int, compute_rec !< Compute interpolations. ! deferred public methods - procedure(compute_stencil_rank_1_interface), pass(self), deferred :: compute_stencil_rank_1 !< Compute interpolation values. - procedure(compute_stencil_rank_2_interface), pass(self), deferred :: compute_stencil_rank_2 !< Compute interpolation values. + procedure(compute_int_interface), pass(self), deferred :: compute_int !< Compute interpolations (interpolate). + procedure(compute_rec_interface), pass(self), deferred :: compute_rec !< Compute interpolations (reconstruct). endtype interpolations_object abstract interface !< Abstract interfaces of [[interpolations_object]]. - pure subroutine compute_stencil_rank_1_interface(self, stencil) - !< Compute interpolations. + pure subroutine compute_int_interface(self, stencil, values) + !< Compute interpolations (interpolate). import :: interpolations_object, RPP - class(interpolations_object), intent(inout) :: self !< Interpolations. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - endsubroutine compute_stencil_rank_1_interface + class(interpolations_object), intent(in) :: self !< Interpolations. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Interpolations values. + endsubroutine compute_int_interface - pure subroutine compute_stencil_rank_2_interface(self, stencil, values) - !< Compute interpolations. + pure subroutine compute_rec_interface(self, stencil, values) + !< Compute interpolations (reconstruct). import :: interpolations_object, RPP class(interpolations_object), intent(in) :: self !< Interpolations. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:, 0:) !< Interpolations values. - endsubroutine compute_stencil_rank_2_interface + endsubroutine compute_rec_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 index d40cbc7..e4f5943 100644 --- a/src/lib/abstract_objects/wenoof_interpolator_object.F90 +++ b/src/lib/abstract_objects/wenoof_interpolator_object.F90 @@ -7,9 +7,9 @@ module wenoof_interpolator_object #else use penf, only: RPP=>R8P #endif -use wenoof_base_object -use wenoof_interpolations_object -use wenoof_weights_object +use wenoof_base_object, only : base_object, base_object_constructor +use wenoof_interpolations_object, only : interpolations_object, interpolations_object_constructor +use wenoof_weights_object, only : weights_object, weights_object_constructor implicit none private @@ -32,52 +32,51 @@ module wenoof_interpolator_object class(weights_object), allocatable :: weights !< Weights of interpolations. contains ! public methods - generic :: interpolate => interpolate_stencil_rank_1_debug, interpolate_stencil_rank_2_debug, & - interpolate_stencil_rank_1_standard, interpolate_stencil_rank_2_standard + generic :: interpolate => interpolate_int_debug, interpolate_rec_debug, & + interpolate_int_standard, interpolate_rec_standard !< Interpolate. ! public deferred methods - procedure(interpolate_stencil_rank_1_debug_interface), pass(self), deferred :: interpolate_stencil_rank_1_debug - procedure(interpolate_stencil_rank_2_debug_interface), pass(self), deferred :: interpolate_stencil_rank_2_debug - procedure(interpolate_stencil_rank_1_standard_interface), pass(self), deferred :: interpolate_stencil_rank_1_standard - procedure(interpolate_stencil_rank_2_standard_interface), pass(self), deferred :: interpolate_stencil_rank_2_standard + procedure(interpolate_int_debug_interface), pass(self), deferred :: interpolate_int_debug !< Interpolate (int) debug. + procedure(interpolate_int_standard_interface), pass(self), deferred :: interpolate_int_standard !< Interpolate (int) standard. + procedure(interpolate_rec_debug_interface), pass(self), deferred :: interpolate_rec_debug !< Interpolate (rec) debug. + procedure(interpolate_rec_standard_interface), pass(self), deferred :: interpolate_rec_standard !< Interpolate (rec) standard. endtype interpolator_object abstract interface !< Abstract interfaces of [[interpolator_object]]. - pure subroutine interpolate_stencil_rank_1_debug_interface(self, stencil, interpolation, si, weights) - !< Interpolate values (providing also debug values). + pure subroutine interpolate_int_debug_interface(self, stencil, interpolation, si, weights) + !< Interpolate (interpolate) values (providing also debug values). import :: interpolator_object, RPP - class(interpolator_object), intent(inout) :: self !< Interpolator. - real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. - real(RPP), intent(out) :: interpolation !< Result of the interpolation. - real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [0:S-1]. - real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [0:S-1]. - endsubroutine interpolate_stencil_rank_1_debug_interface + class(interpolator_object), intent(in) :: self !< Interpolator. + real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. + real(RPP), intent(out) :: interpolation !< Result of the interpolation. + real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [0:S-1]. + real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [0:S-1]. + endsubroutine interpolate_int_debug_interface - pure subroutine interpolate_stencil_rank_2_debug_interface(self, stencil, interpolation, si, weights) - !< Interpolate values (providing also debug values). + pure subroutine interpolate_int_standard_interface(self, stencil, interpolation) + !< Interpolate (interpolate) values (without providing debug values). import :: interpolator_object, RPP - class(interpolator_object), intent(inout) :: self !< Interpolator. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. - real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - endsubroutine interpolate_stencil_rank_2_debug_interface + class(interpolator_object), intent(in) :: self !< Interpolator. + real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. + real(RPP), intent(out) :: interpolation !< Result of the interpolation. + endsubroutine interpolate_int_standard_interface - pure subroutine interpolate_stencil_rank_1_standard_interface(self, stencil, interpolation) - !< Interpolate values (without providing debug values). + pure subroutine interpolate_rec_debug_interface(self, stencil, interpolation, si, weights) + !< Interpolate (reconstruct) values (providing also debug values). import :: interpolator_object, RPP - class(interpolator_object), intent(inout) :: self !< Interpolator. - real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. - real(RPP), intent(out) :: interpolation !< Result of the interpolation. - endsubroutine interpolate_stencil_rank_1_standard_interface + class(interpolator_object), intent(in) :: self !< Interpolator. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. + endsubroutine interpolate_rec_debug_interface - pure subroutine interpolate_stencil_rank_2_standard_interface(self, stencil, interpolation) - !< Interpolate values (without providing debug values). + pure subroutine interpolate_rec_standard_interface(self, stencil, interpolation) + !< Interpolate (reconstruct) values (without providing debug values). import :: interpolator_object, RPP - class(interpolator_object), intent(inout) :: self !< Interpolator. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - endsubroutine interpolate_stencil_rank_2_standard_interface + class(interpolator_object), intent(in) :: self !< Interpolator. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + endsubroutine interpolate_rec_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 index b395b02..edb3d5d 100644 --- a/src/lib/abstract_objects/wenoof_kappa_object.F90 +++ b/src/lib/abstract_objects/wenoof_kappa_object.F90 @@ -7,7 +7,7 @@ module wenoof_kappa_object #else use penf, only: I_P, RPP=>R8P #endif -use wenoof_base_object +use wenoof_base_object, only : base_object, base_object_constructor implicit none private @@ -16,37 +16,34 @@ module wenoof_kappa_object type, extends(base_object_constructor), abstract :: kappa_object_constructor !< Abstract kappa object constructor. - real(RPP), allocatable :: stencil(:) !< Stencil used for interpolation, [1-S:S-1]. - real(RPP) :: x_target !< Coordinate of the interpolation point. endtype kappa_object_constructor type, extends(base_object), abstract :: kappa_object !< Kappa (optimal, linear weights of stencil interpolations) object. - real(RPP), allocatable :: values_rank_1(:) !< Kappa coefficients values [0:S-1]. - real(RPP), allocatable :: values_rank_2(:,:) !< Kappa coefficients values [1:2,0:S-1]. contains ! public methods - generic :: compute => compute_kappa_rec, compute_kappa_int + generic :: compute => compute_rec, compute_int !< Compute kappa. ! deferred public methods - procedure(compute_kappa_rec_interface), pass(self), deferred :: compute_kappa_rec!< Compute interp. - procedure(compute_kappa_int_interface), pass(self), deferred :: compute_kappa_int!< Compute interp. + procedure(compute_int_interface), pass(self), deferred :: compute_int !< Compute kappa (interpolate). + procedure(compute_rec_interface), pass(self), deferred :: compute_rec !< Compute kappa (reconstruct). endtype kappa_object abstract interface !< Abstract interfaces of [[kappa_object]]. - pure subroutine compute_kappa_rec_interface(self) - !< Compute kappa. - import :: kappa_object - class(kappa_object), intent(inout) :: self !< Kappa. - endsubroutine compute_kappa_rec_interface - - pure subroutine compute_kappa_int_interface(self, stencil, x_target) - !< Compute kappa. + pure subroutine compute_int_interface(self, stencil, x_target, values) + !< Compute kappa (interpolate). import :: kappa_object, I_P, RPP - class(kappa_object), intent(inout) :: self !< Kappa. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for interpolation, [1-S:S-1]. - real(RPP), intent(in) :: x_target !< Coordinate of the interpolation point. - endsubroutine compute_kappa_int_interface -endinterface + class(kappa_object), intent(in) :: self !< Kappa. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for interpolation, [1-S:S-1]. + real(RPP), intent(in) :: x_target !< Coordinate of the interpolation point. + real(RPP), intent(out) :: values(0:) !< Kappa values. + endsubroutine compute_int_interface + pure subroutine compute_rec_interface(self, values) + !< Compute kappa (reconstruct). + import :: kappa_object, RPP + class(kappa_object), intent(in) :: self !< Kappa. + real(RPP), intent(out) :: values(1:,0:) !< Kappa values. + endsubroutine compute_rec_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 index 973656a..8eb0442 100644 --- a/src/lib/abstract_objects/wenoof_weights_object.F90 +++ b/src/lib/abstract_objects/wenoof_weights_object.F90 @@ -7,7 +7,7 @@ module wenoof_weights_object #else use penf, only: RPP=>R8P #endif -use wenoof_base_object +use wenoof_base_object, only : base_object, base_object_constructor implicit none private @@ -20,47 +20,47 @@ module wenoof_weights_object type, extends(base_object), abstract :: weights_object !< Weights of stencil interpolations object. - real(RPP), allocatable :: values_rank_1(:) !< Weights values of stencil interpolations [0:S-1]. contains ! public methods - generic :: compute => compute_stencil_rank_1, compute_stencil_rank_2 !< Compute weights. - generic :: smoothness_indicators => smoothness_indicators_of_rank_1, smoothness_indicators_of_rank_2 !< Return IS. + generic :: compute => compute_int, compute_rec !< Compute weights. + generic :: smoothness_indicators => smoothness_indicators_int, smoothness_indicators_rec !< Return IS. ! deferred public methods - procedure(compute_stencil_rank_1_interface), pass(self), deferred :: compute_stencil_rank_1 !< Compute weights. - procedure(compute_stencil_rank_2_interface), pass(self), deferred :: compute_stencil_rank_2 !< Compute weights. - procedure(smoothness_indicators_of_rank_1_interface), pass(self), deferred :: smoothness_indicators_of_rank_1 !< Return IS. - procedure(smoothness_indicators_of_rank_2_interface), pass(self), deferred :: smoothness_indicators_of_rank_2 !< Return IS. + procedure(compute_int_interface), pass(self), deferred :: compute_int !< Compute weights (interp.). + procedure(compute_rec_interface), pass(self), deferred :: compute_rec !< Compute weights (recons.). + procedure(smoothness_indicators_int_interface), pass(self), deferred :: smoothness_indicators_int !< Return IS. + procedure(smoothness_indicators_rec_interface), pass(self), deferred :: smoothness_indicators_rec !< Return IS. endtype weights_object abstract interface !< Abstract interfaces of [[weights_object]]. - pure subroutine compute_stencil_rank_1_interface(self, stencil) - !< Compute beta. + pure subroutine compute_int_interface(self, stencil, values) + !< Compute weights (interpolate). import :: weights_object, RPP - class(weights_object), intent(inout) :: self !< Weights. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - endsubroutine compute_stencil_rank_1_interface + class(weights_object), intent(in) :: self !< Weights. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Weights values of stencil interpolations. + endsubroutine compute_int_interface - pure subroutine compute_stencil_rank_2_interface(self, stencil, values) - !< Compute beta. + pure subroutine compute_rec_interface(self, stencil, values) + !< Compute beta (reconstruct). import :: weights_object, RPP class(weights_object), intent(in) :: self !< Weights. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:,0:) !< Weights values of stencil interpolations. - endsubroutine compute_stencil_rank_2_interface + endsubroutine compute_rec_interface - pure subroutine smoothness_indicators_of_rank_1_interface(self, si) - !< Return smoothness indicators. + pure subroutine smoothness_indicators_int_interface(self, si) + !< Return smoothness indicators (interpolate). import :: weights_object, RPP class(weights_object), intent(in) :: self !< Weights. real(RPP), intent(out) :: si(:) !< Smoothness indicators. - endsubroutine smoothness_indicators_of_rank_1_interface + endsubroutine smoothness_indicators_int_interface - pure subroutine smoothness_indicators_of_rank_2_interface(self, si) - !< Return smoothness indicators. + pure subroutine smoothness_indicators_rec_interface(self, si) + !< Return smoothness indicators (reconstruct). import :: weights_object, RPP - class(weights_object), intent(in) :: self !< Weights. - real(RPP), intent(out) :: si(:,:) !< Smoothness indicators. - endsubroutine smoothness_indicators_of_rank_2_interface + class(weights_object), intent(in) :: self !< Weights. + real(RPP), intent(out) :: si(1:,0:) !< Smoothness indicators. + endsubroutine smoothness_indicators_rec_interface endinterface endmodule wenoof_weights_object diff --git a/src/lib/concrete_objects/wenoof_alpha_int_js.F90 b/src/lib/concrete_objects/wenoof_alpha_int_js.F90 index be87620..8e68113 100644 --- a/src/lib/concrete_objects/wenoof_alpha_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_alpha_int_js.F90 @@ -10,10 +10,8 @@ module wenoof_alpha_int_js #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_alpha_object -use wenoof_base_object -use wenoof_beta_object -use wenoof_kappa_object +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_base_object, only : base_object_constructor implicit none private @@ -31,11 +29,11 @@ module wenoof_alpha_int_js !< 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_interpolate !< Compute alpha (interpolate). - procedure, pass(self) :: compute_reconstruct !< Compute alpha (reconstruct). - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute_int !< Compute alpha (interpolate). + procedure, pass(self) :: compute_rec !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy alpha. endtype alpha_int_js contains @@ -47,15 +45,10 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - associate(val => self%values_rank_1, val_sum => self%values_sum_rank_1) - val = 0._RPP - val_sum = 0._RPP - endassociate endsubroutine create - pure subroutine compute_interpolate(self, beta, kappa, values) - !< Compute alpha. + pure subroutine compute_int(self, beta, kappa, values) + !< Compute alpha (interpolate). class(alpha_int_js), intent(in) :: self !< Alpha coefficient. real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. @@ -65,27 +58,29 @@ pure subroutine compute_interpolate(self, beta, kappa, values) do s1=0, self%S - 1 ! stencil loops values(s1) = kappa(s1) / (self%eps + beta(s1)) ** self%S enddo - endsubroutine compute_interpolate + endsubroutine compute_int - pure subroutine compute_reconstruct(self, beta, kappa, values) - !< Compute alpha. + pure subroutine compute_rec(self, beta, kappa, values) + !< Compute alpha (reconstruct). class(alpha_int_js), intent(in) :: self !< Alpha coefficient. real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + ! empty procedure + endsubroutine compute_rec - ! Empty procedure. - endsubroutine compute_reconstruct + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(alpha_int_js), intent(in) :: self !< Alpha coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. - pure function description(self) result(string) - !< Return alpha string-descripition. - class(alpha_int_js), intent(in) :: self !< Alpha coefficient. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Jiang-Shu alpha coefficients for reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//' - eps = '//trim(str(self%eps)) + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu alpha coefficients object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//' - eps = '//trim(str(self%eps)) endfunction description elemental subroutine destroy(self) @@ -93,6 +88,5 @@ elemental subroutine destroy(self) class(alpha_int_js), intent(inout) :: self !< Alpha. call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) endsubroutine destroy endmodule wenoof_alpha_int_js diff --git a/src/lib/concrete_objects/wenoof_alpha_int_m.F90 b/src/lib/concrete_objects/wenoof_alpha_int_m.F90 new file mode 100644 index 0000000..d471c69 --- /dev/null +++ b/src/lib/concrete_objects/wenoof_alpha_int_m.F90 @@ -0,0 +1,132 @@ +!< Henrick alpha (non linear weights) object. +module wenoof_alpha_int_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 + +#ifdef r16p +use penf, only: I_P, RPP=>R16P, str +#else +use penf, only: I_P, RPP=>R8P, str +#endif +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_alpha_rec_js, only : alpha_rec_js, alpha_rec_js_constructor +use wenoof_alpha_rec_z, only : alpha_rec_z, alpha_rec_z_constructor +use wenoof_base_object, only : base_object_constructor + +implicit none +private +public :: alpha_int_m +public :: alpha_int_m_constructor + +type, extends(alpha_object_constructor) :: alpha_int_m_constructor + !< Henrick alpha (non linear weights) object constructor. + character(len=:), allocatable :: base_type !< Base alpha coefficient type. +endtype alpha_int_m_constructor + +type, extends(alpha_object) :: alpha_int_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_int !< Compute alpha (interpolate). + procedure, pass(self) :: compute_rec !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy alpha. +endtype alpha_int_m + +contains + ! deferred public methods + subroutine create(self, constructor) + !< Create alpha. + class(alpha_int_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_int_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_int(self, beta, kappa, values) + !< Compute alpha (interpolate). + class(alpha_int_m), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + real(RPP) :: alpha_base(0:self%S-1) !< Alpha base coefficients. + real(RPP) :: alpha_base_sum !< Sum of alpha base coefficients. + real(RPP) :: kappa_base !< Kappa base coefficient. + integer(I_P) :: s1 !< Counter. + + call self%alpha_base%compute(beta=beta, kappa=kappa, values=alpha_base) + alpha_base_sum = sum(alpha_base) + do s1=0, self%S - 1 ! stencil loops + kappa_base = alpha_base(s1) / alpha_base_sum + values(s1) = & + (kappa_base * (kappa(s1) + kappa(s1) * kappa(s1) - 3._RPP * kappa(s1) * kappa_base + kappa_base * kappa_base)) / & + (kappa(s1) * kappa(s1) + kappa_base * (1._RPP - 2._RPP * kappa(s1))) + enddo + endsubroutine compute_int + + pure subroutine compute_rec(self, beta, kappa, values) + !< Compute alpha (reconstruct). + class(alpha_int_m), intent(in) :: self !< Alpha coefficient. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + ! empty procedure + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(alpha_int_m), intent(in) :: self !< Alpha coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu alpha coefficients object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//' - eps = '//trim(str(self%eps)) + associate(alpha_base=>self%alpha_base) + select type(alpha_base) + type is(alpha_rec_js) + string = prefix_//string//' - base-mapped-alpha type = Jiang-Shu' + type is(alpha_rec_z) + string = prefix_//string//' - base-mapped-alpha type = Bogeg' + endselect + endassociate + endfunction description + + elemental subroutine destroy(self) + !< Destroy alpha. + class(alpha_int_m), intent(inout) :: self !< Alpha. + + call self%destroy_ + if (allocated(self%alpha_base)) deallocate(self%alpha_base) + endsubroutine destroy +endmodule wenoof_alpha_int_m diff --git a/quarantena/wenoof_alpha_int_z.F90 b/src/lib/concrete_objects/wenoof_alpha_int_z.F90 similarity index 56% rename from quarantena/wenoof_alpha_int_z.F90 rename to src/lib/concrete_objects/wenoof_alpha_int_z.F90 index c7c113a..287146e 100644 --- a/quarantena/wenoof_alpha_int_z.F90 +++ b/src/lib/concrete_objects/wenoof_alpha_int_z.F90 @@ -11,10 +11,8 @@ module wenoof_alpha_int_z #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_alpha_object -use wenoof_base_object -use wenoof_beta_object -use wenoof_kappa_object +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_base_object, only : base_object_constructor implicit none private @@ -31,14 +29,13 @@ module wenoof_alpha_int_z !< @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. - real(RPP), allocatable :: values(:) !< Alpha coefficients [0:S-1]. - real(RPP) :: values_sum !< Sum of alpha coefficients. contains ! public deferred methods - procedure, pass(self) :: create !< Create alpha. - procedure, pass(self) :: compute => compute_alpha_int !< Compute alpha. - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute_int !< Compute alpha (interpolate). + procedure, pass(self) :: compute_rec !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy alpha. endtype alpha_int_z contains ! public deferred methods @@ -49,41 +46,42 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - associate(val => self%values_rank_1, val_sum => self%values_sum_rank_1) - val = 0._RPP - val_sum = 0._RPP - endassociate endsubroutine create - pure subroutine compute_alpha_int(self, beta, kappa) - !< Compute alpha. - class(alpha_int_z), intent(inout) :: self !< Alpha. - class(beta_object), intent(in) :: beta !< Beta. - class(kappa_object), intent(in) :: kappa !< Kappa. - integer(I_P) :: s1 !< Counter. - - associate(val => self%values_rank_1, val_sum => self%values_sum_rank_1) - val_sum = 0._RPP - do s1=0, self%S - 1 ! stencil loops - val = kappa%values_rank_1(s1) * & - ((1._RPP + (tau(S=self%S, beta=beta%values_rank_1) / & - (self%eps + beta%values_rank_1(s1)))) ** (weno_exp(self%S))) - val_sum = val_sum + val(s1) - enddo - endassociate - endsubroutine compute_alpha_int - - pure function description(self) result(string) - !< Return alpha string-descripition. - class(alpha_int_z), intent(in) :: self !< Alpha coefficients. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Borges alpha coefficients for reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//' - eps = '//trim(str(self%eps)) - + pure subroutine compute_int(self, beta, kappa, values) + !< Compute alpha (interpolate). + class(alpha_int_z), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + integer(I_P) :: s1 !< Counter. + + do s1=0, self%S - 1 ! stencil loops + values = kappa(s1) * ((1._RPP + (tau(S=self%S, beta=beta) / (self%eps + beta(s1)))) ** (weno_exp(self%S))) + enddo + endsubroutine compute_int + + pure subroutine compute_rec(self, beta, kappa, values) + !< Compute alpha (reconstruct). + class(alpha_int_z), intent(in) :: self !< Alpha coefficient. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + ! empty procedure + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(alpha_int_z), intent(in) :: self !< Alpha coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu alpha coefficients object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//' - eps = '//trim(str(self%eps)) endfunction description elemental subroutine destroy(self) @@ -91,7 +89,6 @@ elemental subroutine destroy(self) class(alpha_int_z), intent(inout) :: self !< Alpha. call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) endsubroutine destroy ! private non TBP diff --git a/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 b/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 index a409691..9c83b8d 100644 --- a/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_alpha_rec_js.F90 @@ -12,8 +12,6 @@ module wenoof_alpha_rec_js #endif use wenoof_alpha_object, only : alpha_object, alpha_object_constructor use wenoof_base_object, only : base_object_constructor -use wenoof_beta_object, only : beta_object -use wenoof_kappa_object, only : kappa_object implicit none private @@ -31,11 +29,11 @@ module wenoof_alpha_rec_js !< 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_interpolate !< Compute alpha (interpolate). - procedure, pass(self) :: compute_reconstruct !< Compute alpha (reconstruct). - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute_int !< Compute alpha (interpolate). + procedure, pass(self) :: compute_rec !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy alpha. endtype alpha_rec_js contains @@ -49,18 +47,17 @@ subroutine create(self, constructor) call self%create_(constructor=constructor) endsubroutine create - pure subroutine compute_interpolate_interface(self, beta, kappa, values) - !< Compute alpha. + pure subroutine compute_int(self, beta, kappa, values) + !< Compute alpha (interpolate). class(alpha_rec_js), intent(in) :: self !< Alpha coefficient. real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + ! empty procedure + endsubroutine compute_int - ! Empty procedure. - endsubroutine compute_interpolate_interface - - pure subroutine compute_reconstruct(self, beta, kappa, values) - !< Compute alpha. + pure subroutine compute_rec(self, beta, kappa, values) + !< Compute alpha (reconstruct). class(alpha_rec_js), intent(in) :: self !< Alpha coefficient. real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. @@ -72,17 +69,20 @@ pure subroutine compute_reconstruct(self, beta, kappa, values) values(f, s1) = kappa(f, s1) / (self%eps + beta(f, s1)) ** self%S enddo enddo - endsubroutine compute_reconstruct + endsubroutine compute_rec - 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. - character(len=1), parameter :: nl=new_line('a') !< New line char. + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(alpha_rec_js), intent(in) :: self !< Alpha coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. - string = ' Jiang-Shu alpha coefficients for reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//' - eps = '//trim(str(self%eps)) + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu alpha coefficients object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//' - eps = '//trim(str(self%eps)) endfunction description elemental subroutine destroy(self) 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..315346f --- /dev/null +++ b/src/lib/concrete_objects/wenoof_alpha_rec_m.F90 @@ -0,0 +1,135 @@ +!< 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 + +#ifdef r16p +use penf, only: I_P, RPP=>R16P, str +#else +use penf, only: I_P, RPP=>R8P, str +#endif +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_alpha_rec_js, only : alpha_rec_js, alpha_rec_js_constructor +use wenoof_alpha_rec_z, only : alpha_rec_z, alpha_rec_z_constructor +use wenoof_base_object, only : base_object_constructor + +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_int !< Compute alpha (interpolate). + procedure, pass(self) :: compute_rec !< Compute alpha (reconstruct). + 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_int(self, beta, kappa, values) + !< Compute alpha (interpolate). + class(alpha_rec_m), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + ! empty procedure + endsubroutine compute_int + + pure subroutine compute_rec(self, beta, kappa, values) + !< Compute alpha (reconstruct). + class(alpha_rec_m), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + real(RPP) :: alpha_base(1:2,0:self%S-1) !< Alpha base coefficients. + real(RPP) :: alpha_base_sum(1:2) !< Sum of alpha base coefficients. + real(RPP) :: kappa_base !< Kappa base coefficient. + integer(I_P) :: f, s1 !< Counters. + + call self%alpha_base%compute(beta=beta, kappa=kappa, values=alpha_base) + alpha_base_sum(1) = sum(alpha_base(1,:)) + alpha_base_sum(2) = sum(alpha_base(2,:)) + do s1=0, self%S - 1 ! stencil loops + do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + kappa_base = alpha_base(f, s1) / alpha_base_sum(f) + values(f, s1) = & + (kappa_base*(kappa(f, s1) + kappa(f, s1) * kappa(f, s1) - 3._RPP * kappa(f, s1) * kappa_base + kappa_base * kappa_base)) /& + (kappa(f, s1) * kappa(f, s1) + kappa_base * (1._RPP - 2._RPP * kappa(f, s1))) + enddo + enddo + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(alpha_rec_m), intent(in) :: self !< Alpha coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu alpha coefficients object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//' - eps = '//trim(str(self%eps)) + associate(alpha_base=>self%alpha_base) + select type(alpha_base) + type is(alpha_rec_js) + string = prefix_//string//' - base-mapped-alpha type = Jiang-Shu' + type is(alpha_rec_z) + string = prefix_//string//' - base-mapped-alpha type = Bogeg' + endselect + endassociate + 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/quarantena/wenoof_alpha_rec_z.F90 b/src/lib/concrete_objects/wenoof_alpha_rec_z.F90 similarity index 56% rename from quarantena/wenoof_alpha_rec_z.F90 rename to src/lib/concrete_objects/wenoof_alpha_rec_z.F90 index 75ca30d..b2e3fbb 100644 --- a/quarantena/wenoof_alpha_rec_z.F90 +++ b/src/lib/concrete_objects/wenoof_alpha_rec_z.F90 @@ -11,10 +11,8 @@ module wenoof_alpha_rec_z #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_alpha_object -use wenoof_base_object -use wenoof_beta_object -use wenoof_kappa_object +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_base_object, only : base_object_constructor implicit none private @@ -31,14 +29,13 @@ module wenoof_alpha_rec_z !< @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. - real(RPP), allocatable :: values(:,:) !< Alpha coefficients [1:2,0:S-1]. - real(RPP), allocatable :: values_sum(:) !< Sum of alpha coefficients [1:2]. contains ! public deferred methods - procedure, pass(self) :: create !< Create alpha. - procedure, pass(self) :: compute => compute_alpha_rec !< Compute alpha. - procedure, pass(self) :: description !< Return alpha string-description. - procedure, pass(self) :: destroy !< Destroy alpha. + procedure, pass(self) :: create !< Create alpha. + procedure, pass(self) :: compute_int !< Compute alpha (interpolate). + procedure, pass(self) :: compute_rec !< Compute alpha (reconstruct). + procedure, pass(self) :: description !< Return alpha string-description. + procedure, pass(self) :: destroy !< Destroy alpha. endtype alpha_rec_z contains ! public deferred methods @@ -49,44 +46,44 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - allocate(self%values_sum_rank_2(1:2)) - associate(val => self%values_rank_2, val_sum => self%values_sum_rank_2) - val = 0._RPP - val_sum = 0._RPP - endassociate endsubroutine create - pure subroutine compute_alpha_rec(self, beta, kappa) + pure subroutine compute_int(self, beta, kappa, values) + !< Compute alpha (interpolate). + class(alpha_rec_z), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(0:) !< Beta [0:S-1]. + real(RPP), intent(in) :: kappa(0:) !< Kappa [0:S-1]. + real(RPP), intent(out) :: values(0:) !< Alpha values [0:S-1]. + ! empty procedure + endsubroutine compute_int + + pure subroutine compute_rec(self, beta, kappa, values) !< 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. - - associate(val => self%values_rank_2, val_sum => self%values_sum_rank_2) - val_sum = 0._RPP - do s1=0, self%S - 1 ! stencil loops - do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - val(f, s1) = kappa%values_rank_2(f, s1) * & - ((1._RPP + (tau(S=self%S, beta=beta%values_rank_2) / & - (self%eps + beta%values_rank_2(f, s1)))) ** (weno_exp(self%S))) - val_sum(f) = val_sum(f) + val(f, s1) - enddo + class(alpha_rec_z), intent(in) :: self !< Alpha. + real(RPP), intent(in) :: beta(1:,0:) !< Beta [1:2,0:S-1]. + real(RPP), intent(in) :: kappa(1:,0:) !< Kappa [1:2,0:S-1]. + real(RPP), intent(out) :: values(1:,0:) !< Alpha values [1:2,0:S-1]. + integer(I_P) :: f, s1 !< Counters. + + do s1=0, self%S - 1 ! stencil loops + do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + values(f, s1) = kappa(f, s1) * ((1._RPP + (tau(S=self%S, beta=beta) / (self%eps + beta(f, s1)))) ** (weno_exp(self%S))) enddo - endassociate - endsubroutine compute_alpha_rec - - 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. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Borges alpha coefficients for reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//' - eps = '//trim(str(self%eps)) - + enddo + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(alpha_rec_z), intent(in) :: self !< Alpha coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Borges alpha coefficients object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//' - eps = '//trim(str(self%eps)) endfunction description elemental subroutine destroy(self) @@ -94,8 +91,6 @@ elemental subroutine destroy(self) class(alpha_rec_z), intent(inout) :: self !< Alpha. call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) - if (allocated(self%values_sum_rank_2)) deallocate(self%values_sum_rank_2) endsubroutine destroy ! private non TBP diff --git a/src/lib/concrete_objects/wenoof_beta_int_js.F90 b/src/lib/concrete_objects/wenoof_beta_int_js.F90 index 1c1b215..8c9472f 100644 --- a/src/lib/concrete_objects/wenoof_beta_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_beta_int_js.F90 @@ -7,12 +7,12 @@ module wenoof_beta_int_js !< doi:10.1137/070679065. #ifdef r16p -use penf, only: I_P, RPP=>R16P +use penf, only: I_P, RPP=>R16P, str #else -use penf, only: I_P, RPP=>R8P +use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_beta_object +use wenoof_base_object, only : base_object_constructor +use wenoof_beta_object, only : beta_object, beta_object_constructor implicit none private @@ -33,11 +33,11 @@ module wenoof_beta_int_js real(RPP), allocatable :: coef(:,:,:) !< Beta coefficients [0:S-1,0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create beta. - procedure, pass(self) :: compute_stencil_rank_1 !< Compute beta. - procedure, pass(self) :: compute_stencil_rank_2 !< Compute beta. - procedure, pass(self) :: description !< Return beta string-description. - procedure, pass(self) :: destroy !< Destroy beta. + procedure, pass(self) :: create !< Create beta. + procedure, pass(self) :: compute_int !< Compute beta (interpolate). + procedure, pass(self) :: compute_rec !< Compute beta (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy beta. endtype beta_int_js contains @@ -49,8 +49,6 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - self%values_rank_1 = 0._RPP allocate(self%coef(0:self%S - 1, 0:self%S - 1, 0:self%S - 1)) associate(c => self%coef) select case(self%S) @@ -2372,42 +2370,42 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_stencil_rank_1(self, stencil) - !< Compute beta. - class(beta_int_js), intent(inout) :: self !< Beta. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - integer(I_P) :: s1, s2, s3 !< Counters. - - associate(val => self%values_rank_1) - do s1=0, self%S - 1 ! stencils loop - val(s1) = 0._RPP - do s2=0, self%S - 1 - do s3=0, self%S - 1 - val(s1) = val(s1) + self%coef(s3, s2, s1) * stencil(s1 - s3) * stencil(s1 - s2) - enddo + pure subroutine compute_int(self, stencil, values) + !< Compute beta (interpolate). + class(beta_int_js), intent(in) :: self !< Beta. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Beta values [0:S-1]. + integer(I_P) :: s1, s2, s3 !< Counters. + + do s1=0, self%S - 1 ! stencils loop + values(s1) = 0._RPP + do s2=0, self%S - 1 + do s3=0, self%S - 1 + values(s1) = values(s1) + self%coef(s3, s2, s1) * stencil(s1 - s3) * stencil(s1 - s2) enddo enddo - endassociate - endsubroutine compute_stencil_rank_1 + enddo + endsubroutine compute_int - pure subroutine compute_stencil_rank_2(self, stencil, values) - !< Compute beta. + pure subroutine compute_rec(self, stencil, values) + !< Compute beta (reconstruct). class(beta_int_js), intent(in) :: self !< Beta. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:,0:) !< Beta values [1:2,0:S-1]. - - ! Empty subroutine. - endsubroutine compute_stencil_rank_2 - - pure function description(self) result(string) - !< Return beta string-description. - class(beta_int_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_int_js%description to be implemented, do not use!' -#endif + ! empty procedure + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(beta_int_js), intent(in) :: self !< Beta coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu beta coefficients object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) endfunction description elemental subroutine destroy(self) @@ -2415,7 +2413,6 @@ elemental subroutine destroy(self) class(beta_int_js), intent(inout) :: self !< Beta. call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) if (allocated(self%coef)) deallocate(self%coef) endsubroutine destroy endmodule wenoof_beta_int_js diff --git a/src/lib/concrete_objects/wenoof_beta_rec_js.F90 b/src/lib/concrete_objects/wenoof_beta_rec_js.F90 index 2a95068..9299465 100644 --- a/src/lib/concrete_objects/wenoof_beta_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_beta_rec_js.F90 @@ -8,12 +8,12 @@ module wenoof_beta_rec_js !< doi:10.1016/j.jcp.2009.07.039 #ifdef r16p -use penf, only: I_P, RPP=>R16P +use penf, only: I_P, RPP=>R16P, str #else -use penf, only: I_P, RPP=>R8P +use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_beta_object +use wenoof_base_object, only : base_object_constructor +use wenoof_beta_object, only : beta_object, beta_object_constructor implicit none private @@ -35,11 +35,11 @@ module wenoof_beta_rec_js real(RPP), allocatable :: coef(:,:,:) !< Beta coefficients [0:S-1,0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create beta. - procedure, pass(self) :: compute_stencil_rank_1 !< Compute beta. - procedure, pass(self) :: compute_stencil_rank_2 !< Compute beta. - procedure, pass(self) :: description !< Return beta string-description. - procedure, pass(self) :: destroy !< Destroy beta. + procedure, pass(self) :: create !< Create beta. + procedure, pass(self) :: compute_int !< Compute beta (interpolate). + procedure, pass(self) :: compute_rec !< Compute beta (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy beta. endtype beta_rec_js contains @@ -2372,16 +2372,16 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_stencil_rank_1(self, stencil) - !< Compute beta. - class(beta_rec_js), intent(inout) :: self !< Beta. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + pure subroutine compute_int(self, stencil, values) + !< Compute beta (interpolate). + class(beta_rec_js), intent(in) :: self !< Beta. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Beta values [0:S-1]. + ! empty procedure + endsubroutine compute_int - ! Empty routine. - endsubroutine compute_stencil_rank_1 - - pure subroutine compute_stencil_rank_2(self, stencil, values) - !< Compute beta. + pure subroutine compute_rec(self, stencil, values) + !< Compute beta (reconstruct). class(beta_rec_js), intent(in) :: self !< Beta. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:,0:) !< Beta values [1:2,0:S-1]. @@ -2392,22 +2392,24 @@ pure subroutine compute_stencil_rank_2(self, stencil, values) values(f, s1) = 0._RPP do s2=0, self%S - 1 do s3=0, self%S - 1 - values(f, s1) = val(f, s1) + self%coef(s3, s2, s1) * stencil(f, s1 - s3) * stencil(f, s1 - s2) + values(f, s1) = values(f, s1) + self%coef(s3, s2, s1) * stencil(f, s1 - s3) * stencil(f, s1 - s2) enddo enddo enddo enddo - endsubroutine compute_stencil_rank_2 - - 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 + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(beta_rec_js), intent(in) :: self !< Beta coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu beta coefficients object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) endfunction description elemental subroutine destroy(self) diff --git a/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 b/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 index a1b7c02..0e5b1ee 100644 --- a/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_interpolations_int_js.F90 @@ -7,12 +7,12 @@ module wenoof_interpolations_int_js !< doi:10.1137/070679065. #ifdef r16p -use penf, only: I_P, RPP=>R16P +use penf, only : I_P, RPP=>R16P, str #else -use penf, only: I_P, RPP=>R8P +use penf, only : I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_interpolations_object +use wenoof_base_object, only : base_object_constructor +use wenoof_interpolations_object, only : interpolations_object, interpolations_object_constructor implicit none private @@ -29,14 +29,14 @@ module wenoof_interpolations_int_js !< @note The provided interpolations implement the Lagrange interpolations defined in *High Order Weighted Essentially !< Nonoscillatory Schemes for Convection Dominated Problems*, Chi-Wang Shu, SIAM Review, 2009, vol. 51, pp. 82--126, !< doi:10.1137/070679065. - real(RPP), allocatable :: coef(:,:) !< Polynomial coefficients [0:S-1,0:S-1]. + real(RPP), allocatable :: coef(:,:) !< Polynomial coefficients [0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create interpolations. - procedure, pass(self) :: compute_stencil_rank_1 !< Compute interpolations. - procedure, pass(self) :: compute_stencil_rank_2 !< Compute interpolations. - procedure, pass(self) :: description !< Return interpolations string-description. - procedure, pass(self) :: destroy !< Destroy interpolations. + procedure, pass(self) :: create !< Create interpolations. + procedure, pass(self) :: compute_int !< Compute interpolations (interpolate). + procedure, pass(self) :: compute_rec !< Compute interpolations (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy interpolations. endtype interpolations_int_js contains @@ -52,8 +52,6 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - self%values_rank_1 = 0._RPP allocate(self%coef(0:self%S - 1, 0:self%S - 1)) select type(constructor) type is(interpolations_int_js_constructor) @@ -346,41 +344,41 @@ subroutine create(self, constructor) endselect endsubroutine create - pure subroutine compute_stencil_rank_1(self, stencil) - !< Compute interpolations. - class(interpolations_int_js), intent(inout) :: self !< Interpolations. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - integer(I_P) :: s1 !< Counter. - integer(I_P) :: s2 !< Counter. + pure subroutine compute_int(self, stencil, values) + !< Compute interpolations (interpolation). + class(interpolations_int_js), intent(in) :: self !< Interpolations. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Interpolations values. + integer(I_P) :: s1 !< Counter. + integer(I_P) :: s2 !< Counter. - associate(val => self%values_rank_1) - val = 0._RPP + values = 0._RPP do s1=0, self%S - 1 ! stencils loop do s2=0, self%S - 1 ! values loop - val(s1) = val(s1) + self%coef(s2, s1) * stencil(-s2 + s1) + values(s1) = values(s1) + self%coef(s2, s1) * stencil(-s2 + s1) enddo enddo - endassociate - endsubroutine compute_stencil_rank_1 + endsubroutine compute_int - pure subroutine compute_stencil_rank_2(self, stencil, values) - !< Compute interpolations. + pure subroutine compute_rec(self, stencil, values) + !< Compute interpolations (reconstruct). class(interpolations_int_js), intent(in) :: self !< Interpolations. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:, 0:) !< Interpolations values. + ! empty procedure + endsubroutine compute_rec - ! Empty Subroutine. - endsubroutine compute_stencil_rank_2 + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(interpolations_int_js), intent(in) :: self !< Interpolations. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. - pure function description(self) result(string) - !< Return interpolations string-description. - class(interpolations_int_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_int_js%description to be implemented, do not use!' -#endif + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu beta interpolations object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) endfunction description elemental subroutine destroy(self) @@ -388,7 +386,6 @@ elemental subroutine destroy(self) class(interpolations_int_js), intent(inout) :: self !< Interpolations. call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) if (allocated(self%coef)) deallocate(self%coef) endsubroutine destroy endmodule wenoof_interpolations_int_js diff --git a/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 b/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 index b4fcd44..08e6d9d 100644 --- a/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_interpolations_rec_js.F90 @@ -8,12 +8,12 @@ module wenoof_interpolations_rec_js !< doi:10.1016/j.jcp.2009.07.039 #ifdef r16p -use penf, only: I_P, RPP=>R16P +use penf, only: I_P, RPP=>R16P, str #else -use penf, only: I_P, RPP=>R8P +use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_interpolations_object +use wenoof_base_object, only : base_object_constructor +use wenoof_interpolations_object, only : interpolations_object, interpolations_object_constructor implicit none private @@ -35,11 +35,11 @@ module wenoof_interpolations_rec_js real(RPP), allocatable :: coef(:,:,:) !< Polynomial coefficients [1:2,0:S-1,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create interpolations. - procedure, pass(self) :: compute_stencil_rank_1 !< Compute interpolations. - procedure, pass(self) :: compute_stencil_rank_2 !< Compute interpolations. - procedure, pass(self) :: description !< Return interpolations string-description. - procedure, pass(self) :: destroy !< Destroy interpolations. + procedure, pass(self) :: create !< Create interpolations. + procedure, pass(self) :: compute_int !< Compute interpolations (interpolate). + procedure, pass(self) :: compute_rec !< Compute interpolations (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy interpolations. endtype interpolations_rec_js contains @@ -449,16 +449,16 @@ subroutine create(self, constructor) endassociate endsubroutine create - pure subroutine compute_stencil_rank_1(self, stencil) - !< Compute interpolations. - class(interpolations_rec_js), intent(inout) :: self !< Interpolations. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + pure subroutine compute_int(self, stencil, values) + !< Compute interpolations (interpolate). + class(interpolations_rec_js), intent(in) :: self !< Interpolations. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Interpolations values. + ! empty procedure + endsubroutine compute_int - ! Empty Subroutine. - endsubroutine compute_stencil_rank_1 - - pure subroutine compute_stencil_rank_2(self, stencil, values) - !< Compute interpolations. + pure subroutine compute_rec(self, stencil, values) + !< Compute interpolations (reconstruct). class(interpolations_rec_js), intent(in) :: self !< Interpolations. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:, 0:) !< Interpolations values. @@ -474,17 +474,19 @@ pure subroutine compute_stencil_rank_2(self, stencil, values) enddo enddo enddo - endsubroutine compute_stencil_rank_2 + endsubroutine compute_rec - pure function description(self) result(string) - !< Return interpolations string-description. - class(interpolations_rec_js), intent(in) :: self !< Interpolations. - character(len=:), allocatable :: string !< String-description. + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(interpolations_rec_js), intent(in) :: self !< Interpolations. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. -#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 + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu beta interpolations object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) endfunction description elemental subroutine destroy(self) diff --git a/src/lib/concrete_objects/wenoof_interpolator_js.F90 b/src/lib/concrete_objects/wenoof_interpolator_js.F90 index 7aac609..66d64bb 100644 --- a/src/lib/concrete_objects/wenoof_interpolator_js.F90 +++ b/src/lib/concrete_objects/wenoof_interpolator_js.F90 @@ -8,12 +8,12 @@ module wenoof_interpolator_js #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_interpolations_factory -use wenoof_interpolations_object -use wenoof_interpolator_object -use wenoof_weights_factory -use wenoof_weights_object +use wenoof_base_object, only : base_object_constructor +use wenoof_interpolations_factory, only : interpolations_factory +use wenoof_interpolations_object, only : interpolations_object +use wenoof_interpolator_object, only : interpolator_object, interpolator_object_constructor +use wenoof_weights_factory, only : weights_factory +use wenoof_weights_object, only : weights_object implicit none private @@ -34,13 +34,13 @@ module wenoof_interpolator_js !< 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 interpolator. - procedure, pass(self) :: description !< Return interpolator string-description. - procedure, pass(self) :: destroy !< Destroy interpolator. - procedure, pass(self) :: interpolate_stencil_rank_1_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_stencil_rank_2_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_stencil_rank_1_debug !< Interpolate values (providing also debug values). - procedure, pass(self) :: interpolate_stencil_rank_2_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: create !< Create interpolator. + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy interpolator. + procedure, pass(self) :: interpolate_int_debug !< Interpolate values (providing also debug values, interpolate). + procedure, pass(self) :: interpolate_int_standard !< Interpolate values (without providing debug values, interpolate). + procedure, pass(self) :: interpolate_rec_debug !< Interpolate values (providing also debug values, reconstruct). + procedure, pass(self) :: interpolate_rec_standard !< Interpolate values (without providing debug values, reconstruct). endtype interpolator_js contains @@ -61,15 +61,18 @@ subroutine create(self, constructor) endselect endsubroutine create - 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. - - string = 'Jiang-Shu reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//self%weights%description() + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(interpolator_js), intent(in) :: self !< Interpolator. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line character. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu interpolator:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//self%weights%description(prefix=prefix_//' ') endfunction description elemental subroutine destroy(self) @@ -81,51 +84,57 @@ elemental subroutine destroy(self) if (allocated(self%weights)) deallocate(self%weights) endsubroutine destroy - pure subroutine interpolate_stencil_rank_1_debug(self, stencil, interpolation, si, weights) - !< Interpolate values (providing also debug values). - class(interpolator_js), intent(inout) :: self !< Interpolator. - real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation !< Result of the interpolation. - real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. - real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [1:2, 0:S-1]. - - call self%interpolate(stencil=stencil, interpolation=interpolation) - call self%weights%smoothness_indicators_of_rank_1(si=si) - weights = self%weights%values_rank_1 - endsubroutine interpolate_stencil_rank_1_debug - - pure subroutine interpolate_stencil_rank_2_debug(self, stencil, interpolation, si, weights) - !< Interpolate values (providing also debug values). - class(interpolator_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. - real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - - ! Empty subroutine. - endsubroutine interpolate_stencil_rank_2_debug - - pure subroutine interpolate_stencil_rank_1_standard(self, stencil, interpolation) - !< Interpolate values (without providing debug values). - class(interpolator_js), intent(inout) :: self !< Interpolator. - real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation !< Result of the interpolation. - integer(I_P) :: s !< Counters. - - call self%interpolations%compute(stencil=stencil) - call self%weights%compute(stencil=stencil) + pure subroutine interpolate_int_debug(self, stencil, interpolation, si, weights) + !< Interpolate values (providing also debug values, interpolate). + class(interpolator_js), intent(in) :: self !< Interpolator. + real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation !< Result of the interpolation. + real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [1:2, 0:S-1]. + real(RPP) :: interpolations(0:self%S - 1) !< Stencils interpolations. + integer(I_P) :: s !< Counters. + + call self%interpolations%compute(stencil=stencil, values=interpolations) + call self%weights%compute(stencil=stencil, values=weights) + ! call self%weights%smoothness_indicators_of_rank_1(si=si) interpolation = 0._RPP do s=0, self%S - 1 ! stencils loop - interpolation = interpolation + self%weights%values_rank_1(s) * self%interpolations%values_rank_1(s) + interpolation = interpolation + weights(s) * interpolations(s) enddo - endsubroutine interpolate_stencil_rank_1_standard - - pure subroutine interpolate_stencil_rank_2_standard(self, stencil, interpolation) - !< Interpolate values (without providing debug values). - class(interpolator_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - - ! Empty subroutine. - endsubroutine interpolate_stencil_rank_2_standard + endsubroutine interpolate_int_debug + + pure subroutine interpolate_int_standard(self, stencil, interpolation) + !< Interpolate values (without providing debug values, interpolate). + class(interpolator_js), intent(in) :: self !< Interpolator. + real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation !< Result of the interpolation. + real(RPP) :: interpolations(0:self%S - 1) !< Stencils interpolations. + real(RPP) :: weights(0:self%S - 1) !< Weights of stencils interpolations. + integer(I_P) :: s !< Counters. + + call self%interpolations%compute(stencil=stencil, values=interpolations) + call self%weights%compute(stencil=stencil, values=weights) + interpolation = 0._RPP + do s=0, self%S - 1 ! stencils loop + interpolation = interpolation + weights(s) * interpolations(s) + enddo + endsubroutine interpolate_int_standard + + pure subroutine interpolate_rec_debug(self, stencil, interpolation, si, weights) + !< Interpolate values (providing also debug values, reconstruct). + class(interpolator_js), intent(in) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. + ! empty procedure + endsubroutine interpolate_rec_debug + + pure subroutine interpolate_rec_standard(self, stencil, interpolation) + !< Interpolate values (without providing debug values, reconstruct). + class(interpolator_js), intent(in) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + ! empty procedure + endsubroutine interpolate_rec_standard endmodule wenoof_interpolator_js diff --git a/src/lib/concrete_objects/wenoof_kappa_int_js.F90 b/src/lib/concrete_objects/wenoof_kappa_int_js.F90 index 4b5705c..4844598 100644 --- a/src/lib/concrete_objects/wenoof_kappa_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_kappa_int_js.F90 @@ -7,15 +7,15 @@ module wenoof_kappa_int_js !< doi:10.1137/070679065. #ifdef r16p -use penf, only: I_P, RPP=>R16P +use penf, only: I_P, RPP=>R16P, str #else -use penf, only: I_P, RPP=>R8P +use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_interpolations_factory -use wenoof_interpolations_object -use wenoof_interpolations_int_js -use wenoof_kappa_object +use wenoof_base_object, only : base_object_constructor +use wenoof_interpolations_factory, only : interpolations_factory +use wenoof_interpolations_object, only : interpolations_object, interpolations_object_constructor +use wenoof_interpolations_int_js, only : interpolations_int_js +use wenoof_kappa_object, only : kappa_object, kappa_object_constructor implicit none private @@ -24,7 +24,9 @@ module wenoof_kappa_int_js type, extends(kappa_object_constructor) :: kappa_int_js_constructor !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal kappa object constructor. - class(interpolations_object_constructor), allocatable :: interpolations_constructor !< interpolations coefficients constructor. + class(interpolations_object_constructor), allocatable :: interpolations_constructor !< Interpolations coefficients constructor. + real(RPP), allocatable :: stencil(:) !< Stencil used for interpolation, [1-S:S-1]. + real(RPP) :: x_target !< Coordinate of the interpolation point. endtype kappa_int_js_constructor type, extends(kappa_object):: kappa_int_js @@ -33,14 +35,15 @@ module wenoof_kappa_int_js !< @note The provided WENO kappa implements the linear weights defined in *High Order Weighted Essentially !< Nonoscillatory Schemes for Convection Dominated Problems*, Chi-Wang Shu, SIAM Review, 2009, vol. 51, pp. 82--126, !< doi:10.1137/070679065. - class(interpolations_object), allocatable :: interpolations !< interpolations coefficients. + class(interpolations_object), allocatable :: interpolations !< Interpolations object for getting coefficients. + real(RPP), allocatable :: values(:) !< Kappa coefficients values [0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create kappa. - procedure, pass(self) :: compute_kappa_rec !< Compute kappa. - procedure, pass(self) :: compute_kappa_int !< Compute kappa. - procedure, pass(self) :: description !< Return kappa string-description. - procedure, pass(self) :: destroy !< Destroy kappa. + procedure, pass(self) :: create !< Create kappa. + procedure, pass(self) :: compute_int !< Compute kappa (interpolate). + procedure, pass(self) :: compute_rec !< Compute kappa (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy kappa. endtype kappa_int_js contains @@ -55,157 +58,146 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - self%values_rank_1 = 0._RPP + allocate(self%values(0:self%S - 1)) select type(constructor) type is(kappa_int_js_constructor) - associate(interpolations_constructor=>constructor%interpolations_constructor) - call i_factory%create(constructor=interpolations_constructor, object=self%interpolations) - call self%compute(stencil=constructor%stencil, x_target=constructor%x_target) - endassociate + call i_factory%create(constructor=constructor%interpolations_constructor, object=self%interpolations) + call self%compute_int(stencil=constructor%stencil, x_target=constructor%x_target, values=self%values) endselect endsubroutine create - pure subroutine compute_kappa_rec(self) + pure subroutine compute_int(self, stencil, x_target, values) !< Compute kappa. - class(kappa_int_js), intent(inout) :: self !< Kappa. - - ! Empty procedure. - endsubroutine compute_kappa_rec - - pure subroutine compute_kappa_int(self, stencil, x_target) - !< Compute kappa. - class(kappa_int_js), intent(inout) :: self !< Kappa. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for interpolation, [1-S:S-1]. - real(RPP), intent(in) :: x_target !< Coordinate of the interpolation point. - real(RPP) :: coeff(0:2*self%S-2) !< Interpolation coefficients on the whole stencil. - real(RPP) :: coef(0:self%S-1,0:self%S-1) !< Temporary variable. - real(RPP) :: prod !< Temporary variable. - real(RPP) :: coeff_t !< Temporary variable. - real(RPP) :: val_sum !< Temporary variable. - integer(I_P) :: i, j, k !< Counters. + class(kappa_int_js), intent(in) :: self !< Kappa. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for interpolation, [1-S:S-1]. + real(RPP), intent(in) :: x_target !< Coordinate of the interpolation point. + real(RPP), intent(out) :: values(0:) !< Kappa values. + real(RPP) :: coeff(0:2*self%S-2) !< Interpolation coefficients on the whole stencil. + real(RPP) :: coef(0:self%S-1,0:self%S-1) !< Temporary variable. + real(RPP) :: prod !< Temporary variable. + real(RPP) :: coeff_t !< Temporary variable. + real(RPP) :: val_sum !< Temporary variable. + integer(I_P) :: i, j, k !< Counters. - associate(S => self%S, val => self%values_rank_1, interp => self%interpolations) - if(x_target==-0.5_RPP) then - ! left interface (i-1/2) + associate(S=>self%S, interp=>self%interpolations) + if (x_target == -0.5_RPP) then ! left interface (i-1/2) select case(S) case(2) ! 3rd order - val(0) = 3._RPP/4._RPP ! stencil 0 - val(1) = 1._RPP/4._RPP ! stencil 1 + values(0) = 3._RPP/4._RPP ! stencil 0 + values(1) = 1._RPP/4._RPP ! stencil 1 case(3) ! 5th order - val(0) = 5._RPP/16._RPP ! stencil 0 - val(1) = 5._RPP/8._RPP ! stencil 1 - val(2) = 1._RPP/16._RPP ! stencil 2 + values(0) = 5._RPP/16._RPP ! stencil 0 + values(1) = 5._RPP/8._RPP ! stencil 1 + values(2) = 1._RPP/16._RPP ! stencil 2 case(4) ! 7th order - val(0) = 7._RPP/64._RPP ! stencil 0 - val(1) = 35._RPP/64._RPP ! stencil 1 - val(2) = 21._RPP/64._RPP ! stencil 2 - val(3) = 1._RPP/64._RPP ! stencil 3 + values(0) = 7._RPP/64._RPP ! stencil 0 + values(1) = 35._RPP/64._RPP ! stencil 1 + values(2) = 21._RPP/64._RPP ! stencil 2 + values(3) = 1._RPP/64._RPP ! stencil 3 case(5) ! 9th order - val(0) = 9._RPP/256._RPP ! stencil 0 - val(1) = 21._RPP/64._RPP ! stencil 1 - val(2) = 63._RPP/128._RPP ! stencil 2 - val(3) = 9._RPP/64._RPP ! stencil 3 - val(4) = 1._RPP/256._RPP ! stencil 4 + values(0) = 9._RPP/256._RPP ! stencil 0 + values(1) = 21._RPP/64._RPP ! stencil 1 + values(2) = 63._RPP/128._RPP ! stencil 2 + values(3) = 9._RPP/64._RPP ! stencil 3 + values(4) = 1._RPP/256._RPP ! stencil 4 case(6) ! 11th order - val(0) = 11._RPP/1024._RPP ! stencil 0 - val(1) = 165._RPP/1024._RPP ! stencil 1 - val(2) = 231._RPP/512._RPP ! stencil 2 - val(3) = 165._RPP/512._RPP ! stencil 3 - val(4) = 55._RPP/1024._RPP ! stencil 4 - val(5) = 1._RPP/1024._RPP ! stencil 5 + values(0) = 11._RPP/1024._RPP ! stencil 0 + values(1) = 165._RPP/1024._RPP ! stencil 1 + values(2) = 231._RPP/512._RPP ! stencil 2 + values(3) = 165._RPP/512._RPP ! stencil 3 + values(4) = 55._RPP/1024._RPP ! stencil 4 + values(5) = 1._RPP/1024._RPP ! stencil 5 case(7) ! 13th order - val(0) = 13._RPP/4096._RPP ! stencil 0 - val(1) = 143._RPP/2048._RPP ! stencil 1 - val(2) = 1287._RPP/4096._RPP ! stencil 2 - val(3) = 429._RPP/1024._RPP ! stencil 3 - val(4) = 179._RPP/1024._RPP ! stencil 4 - val(5) = 39._RPP/2048._RPP ! stencil 5 - val(6) = 1._RPP/4096._RPP ! stencil 6 + values(0) = 13._RPP/4096._RPP ! stencil 0 + values(1) = 143._RPP/2048._RPP ! stencil 1 + values(2) = 1287._RPP/4096._RPP ! stencil 2 + values(3) = 429._RPP/1024._RPP ! stencil 3 + values(4) = 179._RPP/1024._RPP ! stencil 4 + values(5) = 39._RPP/2048._RPP ! stencil 5 + values(6) = 1._RPP/4096._RPP ! stencil 6 case(8) ! 15th order - val(0) = 15._RPP/16384._RPP ! stencil 0 - val(1) = 455._RPP/16384._RPP ! stencil 1 - val(2) = 3003._RPP/16384._RPP ! stencil 2 - val(3) = 6435._RPP/16384._RPP ! stencil 3 - val(4) = 5005._RPP/16384._RPP ! stencil 4 - val(5) = 1365._RPP/16384._RPP ! stencil 5 - val(6) = 105._RPP/16384._RPP ! stencil 6 - val(7) = 1._RPP/16384._RPP ! stencil 7 + values(0) = 15._RPP/16384._RPP ! stencil 0 + values(1) = 455._RPP/16384._RPP ! stencil 1 + values(2) = 3003._RPP/16384._RPP ! stencil 2 + values(3) = 6435._RPP/16384._RPP ! stencil 3 + values(4) = 5005._RPP/16384._RPP ! stencil 4 + values(5) = 1365._RPP/16384._RPP ! stencil 5 + values(6) = 105._RPP/16384._RPP ! stencil 6 + values(7) = 1._RPP/16384._RPP ! stencil 7 case(9) ! 17th order - val(0) = 17._RPP/65536._RPP ! stencil 0 - val(1) = 85._RPP/8192._RPP ! stencil 1 - val(2) = 1547._RPP/16384._RPP ! stencil 2 - val(3) = 2431._RPP/8192._RPP ! stencil 3 - val(4) = 12155._RPP/32768._RPP ! stencil 4 - val(5) = 1547._RPP/8192._RPP ! stencil 5 - val(6) = 595._RPP/16384._RPP ! stencil 6 - val(7) = 17._RPP/8192._RPP ! stencil 7 - val(8) = 1._RPP/65536._RPP ! stencil 8 + values(0) = 17._RPP/65536._RPP ! stencil 0 + values(1) = 85._RPP/8192._RPP ! stencil 1 + values(2) = 1547._RPP/16384._RPP ! stencil 2 + values(3) = 2431._RPP/8192._RPP ! stencil 3 + values(4) = 12155._RPP/32768._RPP ! stencil 4 + values(5) = 1547._RPP/8192._RPP ! stencil 5 + values(6) = 595._RPP/16384._RPP ! stencil 6 + values(7) = 17._RPP/8192._RPP ! stencil 7 + values(8) = 1._RPP/65536._RPP ! stencil 8 endselect - elseif(x_target==0.5_RPP) then - ! right interface (i+1/2) + elseif(x_target == 0.5_RPP) then ! right interface (i+1/2) select case(S) case(2) ! 3rd order - val(0) = 1._RPP/4._RPP ! stencil 0 - val(1) = 3._RPP/4._RPP ! stencil 1 + values(0) = 1._RPP/4._RPP ! stencil 0 + values(1) = 3._RPP/4._RPP ! stencil 1 case(3) ! 5th order - val(0) = 1._RPP/16._RPP ! stencil 0 - val(1) = 5._RPP/8._RPP ! stencil 1 - val(2) = 5._RPP/16._RPP ! stencil 2 + values(0) = 1._RPP/16._RPP ! stencil 0 + values(1) = 5._RPP/8._RPP ! stencil 1 + values(2) = 5._RPP/16._RPP ! stencil 2 case(4) ! 7th order - val(0) = 1._RPP/64._RPP ! stencil 0 - val(1) = 21._RPP/64._RPP ! stencil 1 - val(2) = 35._RPP/64._RPP ! stencil 2 - val(3) = 7._RPP/64._RPP ! stencil 3 + values(0) = 1._RPP/64._RPP ! stencil 0 + values(1) = 21._RPP/64._RPP ! stencil 1 + values(2) = 35._RPP/64._RPP ! stencil 2 + values(3) = 7._RPP/64._RPP ! stencil 3 case(5) ! 9th order - val(0) = 1._RPP/256._RPP ! stencil 0 - val(1) = 9._RPP/64._RPP ! stencil 1 - val(2) = 63._RPP/128._RPP ! stencil 2 - val(3) = 21._RPP/64._RPP ! stencil 3 - val(4) = 9._RPP/256._RPP ! stencil 4 + values(0) = 1._RPP/256._RPP ! stencil 0 + values(1) = 9._RPP/64._RPP ! stencil 1 + values(2) = 63._RPP/128._RPP ! stencil 2 + values(3) = 21._RPP/64._RPP ! stencil 3 + values(4) = 9._RPP/256._RPP ! stencil 4 case(6) ! 11th order - val(0) = 1._RPP/1024._RPP ! stencil 0 - val(1) = 55._RPP/1024._RPP ! stencil 1 - val(2) = 165._RPP/512._RPP ! stencil 2 - val(3) = 231._RPP/512._RPP ! stencil 3 - val(4) = 165._RPP/1024._RPP ! stencil 4 - val(5) = 11._RPP/1024._RPP ! stencil 5 + values(0) = 1._RPP/1024._RPP ! stencil 0 + values(1) = 55._RPP/1024._RPP ! stencil 1 + values(2) = 165._RPP/512._RPP ! stencil 2 + values(3) = 231._RPP/512._RPP ! stencil 3 + values(4) = 165._RPP/1024._RPP ! stencil 4 + values(5) = 11._RPP/1024._RPP ! stencil 5 case(7) ! 13th order - val(0) = 1._RPP/4096._RPP ! stencil 0 - val(1) = 39._RPP/2048._RPP ! stencil 1 - val(2) = 179._RPP/1024._RPP ! stencil 2 - val(3) = 429._RPP/1024._RPP ! stencil 3 - val(4) = 1287._RPP/4096._RPP ! stencil 4 - val(5) = 143._RPP/2048._RPP ! stencil 5 - val(6) = 13._RPP/4096._RPP ! stencil 6 + values(0) = 1._RPP/4096._RPP ! stencil 0 + values(1) = 39._RPP/2048._RPP ! stencil 1 + values(2) = 179._RPP/1024._RPP ! stencil 2 + values(3) = 429._RPP/1024._RPP ! stencil 3 + values(4) = 1287._RPP/4096._RPP ! stencil 4 + values(5) = 143._RPP/2048._RPP ! stencil 5 + values(6) = 13._RPP/4096._RPP ! stencil 6 case(8) ! 15th order - val(0) = 1._RPP/16384._RPP ! stencil 0 - val(1) = 105._RPP/16384._RPP ! stencil 1 - val(2) = 1365._RPP/16384._RPP ! stencil 2 - val(3) = 5005._RPP/16384._RPP ! stencil 3 - val(4) = 6435._RPP/16384._RPP ! stencil 4 - val(5) = 3003._RPP/16384._RPP ! stencil 5 - val(6) = 455._RPP/16384._RPP ! stencil 6 - val(7) = 15._RPP/16384._RPP ! stencil 7 + values(0) = 1._RPP/16384._RPP ! stencil 0 + values(1) = 105._RPP/16384._RPP ! stencil 1 + values(2) = 1365._RPP/16384._RPP ! stencil 2 + values(3) = 5005._RPP/16384._RPP ! stencil 3 + values(4) = 6435._RPP/16384._RPP ! stencil 4 + values(5) = 3003._RPP/16384._RPP ! stencil 5 + values(6) = 455._RPP/16384._RPP ! stencil 6 + values(7) = 15._RPP/16384._RPP ! stencil 7 case(9) ! 17th order - val(0) = 1._RPP/65536._RPP ! stencil 0 - val(1) = 17._RPP/8192._RPP ! stencil 1 - val(2) = 595._RPP/16384._RPP ! stencil 2 - val(3) = 1547._RPP/8192._RPP ! stencil 3 - val(4) = 12155._RPP/32768._RPP ! stencil 4 - val(5) = 2431._RPP/8192._RPP ! stencil 5 - val(6) = 1547._RPP/16384._RPP ! stencil 6 - val(7) = 85._RPP/8192._RPP ! stencil 7 - val(8) = 17._RPP/65536._RPP ! stencil 8 + values(0) = 1._RPP/65536._RPP ! stencil 0 + values(1) = 17._RPP/8192._RPP ! stencil 1 + values(2) = 595._RPP/16384._RPP ! stencil 2 + values(3) = 1547._RPP/8192._RPP ! stencil 3 + values(4) = 12155._RPP/32768._RPP ! stencil 4 + values(5) = 2431._RPP/8192._RPP ! stencil 5 + values(6) = 1547._RPP/16384._RPP ! stencil 6 + values(7) = 85._RPP/8192._RPP ! stencil 7 + values(8) = 17._RPP/65536._RPP ! stencil 8 endselect - elseif(x_target==0._RPP) then - val = 1._RPP / S + elseif(x_target == 0._RPP) then + values = 1._RPP / S else ! internal point val_sum = 0._RPP - do j=0,2*S-3 !values loop + do j=0, 2 * S - 3 ! values loop prod = 1._RPP - do i=0,2*S-2 + do i=0, 2 * S - 2 if (i==j) cycle prod = prod * ((x_target - stencil(-S+i+1)) / (stencil(-S+j+1) - stencil(-S+i+1))) enddo @@ -218,34 +210,43 @@ pure subroutine compute_kappa_int(self, stencil, x_target) val_sum = 0._RPP do k=0,S-1 do j=0,S-1 - coef(j,k) = interp%coef(S-1-j,S-1-k) + coef(j, k) = interp%coef(S-1-j,S-1-k) enddo enddo do j = 0,S-2 coeff_t = 0._RPP k = j do i = 0,j-1 - coeff_t = coeff_t + val(i) * coef(k,i) + coeff_t = coeff_t + values(i) * coef(k,i) k = k - 1 enddo - val(j) = (coeff(j) - coeff_t) / coef(0,j) - val_sum = val_sum + val(j) + values(j) = (coeff(j) - coeff_t) / coef(0,j) + val_sum = val_sum + values(j) enddo - val(S-1) = 1._RPP - val_sum + values(S-1) = 1._RPP - val_sum endselect endif endassociate - endsubroutine compute_kappa_int + endsubroutine compute_int - pure function description(self) result(string) - !< Return string-description of kappa. - class(kappa_int_js), intent(in) :: self !< Kappa. - character(len=:), allocatable :: string !< String-description. + pure subroutine compute_rec(self, values) + !< Compute kappa (reconstruct). + class(kappa_int_js), intent(in) :: self !< Kappa. + real(RPP), intent(out) :: values(1:,0:) !< Kappa values. + ! empty procedure + endsubroutine compute_rec -#ifndef DEBUG - ! error stop in pure procedure is a F2015 feature not yet supported in debug mode - error stop 'kappa_int_js%description to be implemented, do not use!' -#endif + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(kappa_int_js), intent(in) :: self !< Kappa coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu kappa coefficients object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) endfunction description elemental subroutine destroy(self) @@ -253,6 +254,6 @@ elemental subroutine destroy(self) class(kappa_int_js), intent(inout) :: self !< Kappa. call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) + if (allocated(self%values)) deallocate(self%values) endsubroutine destroy endmodule wenoof_kappa_int_js diff --git a/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 b/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 index fe7a975..e632bd1 100644 --- a/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_kappa_rec_js.F90 @@ -8,12 +8,12 @@ module wenoof_kappa_rec_js !< doi:10.1016/j.jcp.2009.07.039 #ifdef r16p -use penf, only: I_P, RPP=>R16P +use penf, only: I_P, RPP=>R16P, str #else -use penf, only: I_P, RPP=>R8P +use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_kappa_object +use wenoof_base_object, only : base_object_constructor +use wenoof_kappa_object, only : kappa_object, kappa_object_constructor implicit none private @@ -34,11 +34,11 @@ module wenoof_kappa_rec_js real(RPP), allocatable :: values(:,:) !< Kappa coefficients values [1:2,0:S-1]. contains ! public deferred methods - procedure, pass(self) :: create !< Create kappa. - procedure, pass(self) :: compute_kappa_rec !< Compute kappa. - procedure, pass(self) :: compute_kappa_int !< Compute kappa. - procedure, pass(self) :: description !< Return kappa string-description. - procedure, pass(self) :: destroy !< Destroy kappa. + procedure, pass(self) :: create !< Create kappa. + procedure, pass(self) :: compute_int !< Compute kappa (interpolate). + procedure, pass(self) :: compute_rec !< Compute kappa (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy kappa. endtype kappa_rec_js contains @@ -52,151 +52,151 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_2(1:2, 0:self%S - 1)) - self%values_rank_2 = 0._RPP - call self%compute + allocate(self%values(1:2, 0:self%S - 1)) + call self%compute_rec(values=self%values) endsubroutine create - pure subroutine compute_kappa_rec(self) - !< Compute kappa. - class(kappa_rec_js), intent(inout) :: self !< Kappa. - - associate(val => self%values_rank_2) - select case(self%S) - case(2) ! 3rd order - ! 1 => left interface (i-1/2) - val(1, 0) = 2._RPP/3._RPP ! stencil 0 - val(1, 1) = 1._RPP/3._RPP ! stencil 1 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/3._RPP ! stencil 0 - val(2, 1) = 2._RPP/3._RPP ! stencil 1 - case(3) ! 5th order - ! 1 => left interface (i-1/2) - val(1, 0) = 0.3_RPP ! stencil 0 - val(1, 1) = 0.6_RPP ! stencil 1 - val(1, 2) = 0.1_RPP ! stencil 2 - ! 2 => right interface (i+1/2) - val(2, 0) = 0.1_RPP ! stencil 0 - val(2, 1) = 0.6_RPP ! stencil 1 - val(2, 2) = 0.3_RPP ! stencil 2 - case(4) ! 7th order - ! 1 => left interface (i-1/2) - val(1, 0) = 4._RPP/35._RPP ! stencil 0 - val(1, 1) = 18._RPP/35._RPP ! stencil 1 - val(1, 2) = 12._RPP/35._RPP ! stencil 2 - val(1, 3) = 1._RPP/35._RPP ! stencil 3 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/35._RPP ! stencil 0 - val(2, 1) = 12._RPP/35._RPP ! stencil 1 - val(2, 2) = 18._RPP/35._RPP ! stencil 2 - val(2, 3) = 4._RPP/35._RPP ! stencil 3 - case(5) ! 9th order - ! 1 => left interface (i-1/2) - val(1, 0) = 5._RPP/126._RPP ! stencil 0 - val(1, 1) = 20._RPP/63._RPP ! stencil 1 - val(1, 2) = 10._RPP/21._RPP ! stencil 2 - val(1, 3) = 10._RPP/63._RPP ! stencil 3 - val(1, 4) = 1._RPP/126._RPP ! stencil 4 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/126._RPP ! stencil 0 - val(2, 1) = 10._RPP/63._RPP ! stencil 1 - val(2, 2) = 10._RPP/21._RPP ! stencil 2 - val(2, 3) = 20._RPP/63._RPP ! stencil 3 - val(2, 4) = 5._RPP/126._RPP ! stencil 4 - case(6) ! 11th order - ! 1 => left interface (i-1/2) - val(1, 0) = 1._RPP/77._RPP ! stencil 0 - val(1, 1) = 25._RPP/154._RPP ! stencil 1 - val(1, 2) = 100._RPP/231._RPP ! stencil 2 - val(1, 3) = 25._RPP/77._RPP ! stencil 3 - val(1, 4) = 5._RPP/77._RPP ! stencil 4 - val(1, 5) = 1._RPP/462._RPP ! stencil 5 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/462._RPP ! stencil 0 - val(2, 1) = 5._RPP/77._RPP ! stencil 1 - val(2, 2) = 25._RPP/77._RPP ! stencil 2 - val(2, 3) = 100._RPP/231._RPP ! stencil 3 - val(2, 4) = 25._RPP/154._RPP ! stencil 4 - val(2, 5) = 1._RPP/77._RPP ! stencil 5 - case(7) ! 13th order - ! 1 => left interface (i-1/2) - val(1, 0) = 7._RPP/1716._RPP ! stencil 0 - val(1, 1) = 21._RPP/286._RPP ! stencil 1 - val(1, 2) = 175._RPP/572._RPP ! stencil 2 - val(1, 3) = 175._RPP/429._RPP ! stencil 3 - val(1, 4) = 105._RPP/572._RPP ! stencil 4 - val(1, 5) = 7._RPP/286._RPP ! stencil 5 - val(1, 6) = 1._RPP/1716._RPP ! stencil 6 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/1716._RPP ! stencil 0 - val(2, 1) = 7._RPP/286._RPP ! stencil 1 - val(2, 2) = 105._RPP/572._RPP ! stencil 2 - val(2, 3) = 175._RPP/429._RPP ! stencil 3 - val(2, 4) = 175._RPP/572._RPP ! stencil 4 - val(2, 5) = 21._RPP/286._RPP ! stencil 5 - val(2, 6) = 7._RPP/1716._RPP ! stencil 6 - case(8) ! 15th order - ! 1 => left interface (i-1/2) - val(1, 0) = 8._RPP/6435._RPP ! stencil 0 - val(1, 1) = 196._RPP/6435._RPP ! stencil 1 - val(1, 2) = 392._RPP/2145._RPP ! stencil 2 - val(1, 3) = 490._RPP/1287._RPP ! stencil 3 - val(1, 4) = 392._RPP/1287._RPP ! stencil 4 - val(1, 5) = 196._RPP/2145._RPP ! stencil 5 - val(1, 6) = 56._RPP/6435._RPP ! stencil 6 - val(1, 7) = 1._RPP/6435._RPP ! stencil 7 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/6435._RPP ! stencil 0 - val(2, 1) = 56._RPP/6435._RPP ! stencil 1 - val(2, 2) = 196._RPP/2145._RPP ! stencil 2 - val(2, 3) = 392._RPP/1287._RPP ! stencil 3 - val(2, 4) = 490._RPP/1287._RPP ! stencil 4 - val(2, 5) = 392._RPP/2145._RPP ! stencil 5 - val(2, 6) = 196._RPP/6435._RPP ! stencil 6 - val(2, 7) = 8._RPP/6435._RPP ! stencil 7 - case(9) ! 17th order - ! 1 => left interface (i-1/2) - val(1, 0) = 9._RPP/24310._RPP ! stencil 0 - val(1, 1) = 144._RPP/12155._RPP ! stencil 1 - val(1, 2) = 1176._RPP/12155._RPP ! stencil 2 - val(1, 3) = 3528._RPP/12155._RPP ! stencil 3 - val(1, 4) = 882._RPP/2431._RPP ! stencil 4 - val(1, 5) = 2352._RPP/12155._RPP ! stencil 5 - val(1, 6) = 504._RPP/12155._RPP ! stencil 6 - val(1, 7) = 36._RPP/12155._RPP ! stencil 7 - val(1, 8) = 1._RPP/24310._RPP ! stencil 8 - ! 2 => right interface (i+1/2) - val(2, 0) = 1._RPP/24310._RPP ! stencil 0 - val(2, 1) = 36._RPP/12155._RPP ! stencil 1 - val(2, 2) = 504._RPP/12155._RPP ! stencil 2 - val(2, 3) = 2352._RPP/12155._RPP ! stencil 3 - val(2, 4) = 882._RPP/2431._RPP ! stencil 4 - val(2, 5) = 3528._RPP/12155._RPP ! stencil 5 - val(2, 6) = 1176._RPP/12155._RPP ! stencil 6 - val(2, 7) = 144._RPP/12155._RPP ! stencil 7 - val(2, 8) = 9._RPP/24310._RPP ! stencil 8 - endselect - endassociate - endsubroutine compute_kappa_rec + pure subroutine compute_int(self, stencil, x_target, values) + !< Compute kappa (interpolate). + class(kappa_rec_js), intent(in) :: self !< Kappa. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for interpolation, [1-S:S-1]. + real(RPP), intent(in) :: x_target !< Coordinate of the interpolation point. + real(RPP), intent(out) :: values(0:) !< Kappa values. + ! empty procedure + endsubroutine compute_int - pure subroutine compute_kappa_int(self, stencil, x_target) - !< Compute kappa. - class(kappa_rec_js), intent(inout) :: self !< Kappa. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for interpolation, [1-S:S-1]. - real(RPP), intent(in) :: x_target !< Coordinate of the interpolation point. + pure subroutine compute_rec(self, values) + !< Compute kappa (reconstruct). + class(kappa_rec_js), intent(in) :: self !< Kappa. + real(RPP), intent(out) :: values(1:,0:) !< Kappa values. - ! Empty Subroutine - endsubroutine compute_kappa_int + select case(self%S) + case(2) ! 3rd order + ! 1 => left interface (i-1/2) + values(1, 0) = 2._RPP/3._RPP ! stencil 0 + values(1, 1) = 1._RPP/3._RPP ! stencil 1 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/3._RPP ! stencil 0 + values(2, 1) = 2._RPP/3._RPP ! stencil 1 + case(3) ! 5th order + ! 1 => left interface (i-1/2) + values(1, 0) = 0.3_RPP ! stencil 0 + values(1, 1) = 0.6_RPP ! stencil 1 + values(1, 2) = 0.1_RPP ! stencil 2 + ! 2 => right interface (i+1/2) + values(2, 0) = 0.1_RPP ! stencil 0 + values(2, 1) = 0.6_RPP ! stencil 1 + values(2, 2) = 0.3_RPP ! stencil 2 + case(4) ! 7th order + ! 1 => left interface (i-1/2) + values(1, 0) = 4._RPP/35._RPP ! stencil 0 + values(1, 1) = 18._RPP/35._RPP ! stencil 1 + values(1, 2) = 12._RPP/35._RPP ! stencil 2 + values(1, 3) = 1._RPP/35._RPP ! stencil 3 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/35._RPP ! stencil 0 + values(2, 1) = 12._RPP/35._RPP ! stencil 1 + values(2, 2) = 18._RPP/35._RPP ! stencil 2 + values(2, 3) = 4._RPP/35._RPP ! stencil 3 + case(5) ! 9th order + ! 1 => left interface (i-1/2) + values(1, 0) = 5._RPP/126._RPP ! stencil 0 + values(1, 1) = 20._RPP/63._RPP ! stencil 1 + values(1, 2) = 10._RPP/21._RPP ! stencil 2 + values(1, 3) = 10._RPP/63._RPP ! stencil 3 + values(1, 4) = 1._RPP/126._RPP ! stencil 4 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/126._RPP ! stencil 0 + values(2, 1) = 10._RPP/63._RPP ! stencil 1 + values(2, 2) = 10._RPP/21._RPP ! stencil 2 + values(2, 3) = 20._RPP/63._RPP ! stencil 3 + values(2, 4) = 5._RPP/126._RPP ! stencil 4 + case(6) ! 11th order + ! 1 => left interface (i-1/2) + values(1, 0) = 1._RPP/77._RPP ! stencil 0 + values(1, 1) = 25._RPP/154._RPP ! stencil 1 + values(1, 2) = 100._RPP/231._RPP ! stencil 2 + values(1, 3) = 25._RPP/77._RPP ! stencil 3 + values(1, 4) = 5._RPP/77._RPP ! stencil 4 + values(1, 5) = 1._RPP/462._RPP ! stencil 5 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/462._RPP ! stencil 0 + values(2, 1) = 5._RPP/77._RPP ! stencil 1 + values(2, 2) = 25._RPP/77._RPP ! stencil 2 + values(2, 3) = 100._RPP/231._RPP ! stencil 3 + values(2, 4) = 25._RPP/154._RPP ! stencil 4 + values(2, 5) = 1._RPP/77._RPP ! stencil 5 + case(7) ! 13th order + ! 1 => left interface (i-1/2) + values(1, 0) = 7._RPP/1716._RPP ! stencil 0 + values(1, 1) = 21._RPP/286._RPP ! stencil 1 + values(1, 2) = 175._RPP/572._RPP ! stencil 2 + values(1, 3) = 175._RPP/429._RPP ! stencil 3 + values(1, 4) = 105._RPP/572._RPP ! stencil 4 + values(1, 5) = 7._RPP/286._RPP ! stencil 5 + values(1, 6) = 1._RPP/1716._RPP ! stencil 6 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/1716._RPP ! stencil 0 + values(2, 1) = 7._RPP/286._RPP ! stencil 1 + values(2, 2) = 105._RPP/572._RPP ! stencil 2 + values(2, 3) = 175._RPP/429._RPP ! stencil 3 + values(2, 4) = 175._RPP/572._RPP ! stencil 4 + values(2, 5) = 21._RPP/286._RPP ! stencil 5 + values(2, 6) = 7._RPP/1716._RPP ! stencil 6 + case(8) ! 15th order + ! 1 => left interface (i-1/2) + values(1, 0) = 8._RPP/6435._RPP ! stencil 0 + values(1, 1) = 196._RPP/6435._RPP ! stencil 1 + values(1, 2) = 392._RPP/2145._RPP ! stencil 2 + values(1, 3) = 490._RPP/1287._RPP ! stencil 3 + values(1, 4) = 392._RPP/1287._RPP ! stencil 4 + values(1, 5) = 196._RPP/2145._RPP ! stencil 5 + values(1, 6) = 56._RPP/6435._RPP ! stencil 6 + values(1, 7) = 1._RPP/6435._RPP ! stencil 7 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/6435._RPP ! stencil 0 + values(2, 1) = 56._RPP/6435._RPP ! stencil 1 + values(2, 2) = 196._RPP/2145._RPP ! stencil 2 + values(2, 3) = 392._RPP/1287._RPP ! stencil 3 + values(2, 4) = 490._RPP/1287._RPP ! stencil 4 + values(2, 5) = 392._RPP/2145._RPP ! stencil 5 + values(2, 6) = 196._RPP/6435._RPP ! stencil 6 + values(2, 7) = 8._RPP/6435._RPP ! stencil 7 + case(9) ! 17th order + ! 1 => left interface (i-1/2) + values(1, 0) = 9._RPP/24310._RPP ! stencil 0 + values(1, 1) = 144._RPP/12155._RPP ! stencil 1 + values(1, 2) = 1176._RPP/12155._RPP ! stencil 2 + values(1, 3) = 3528._RPP/12155._RPP ! stencil 3 + values(1, 4) = 882._RPP/2431._RPP ! stencil 4 + values(1, 5) = 2352._RPP/12155._RPP ! stencil 5 + values(1, 6) = 504._RPP/12155._RPP ! stencil 6 + values(1, 7) = 36._RPP/12155._RPP ! stencil 7 + values(1, 8) = 1._RPP/24310._RPP ! stencil 8 + ! 2 => right interface (i+1/2) + values(2, 0) = 1._RPP/24310._RPP ! stencil 0 + values(2, 1) = 36._RPP/12155._RPP ! stencil 1 + values(2, 2) = 504._RPP/12155._RPP ! stencil 2 + values(2, 3) = 2352._RPP/12155._RPP ! stencil 3 + values(2, 4) = 882._RPP/2431._RPP ! stencil 4 + values(2, 5) = 3528._RPP/12155._RPP ! stencil 5 + values(2, 6) = 1176._RPP/12155._RPP ! stencil 6 + values(2, 7) = 144._RPP/12155._RPP ! stencil 7 + values(2, 8) = 9._RPP/24310._RPP ! stencil 8 + endselect + endsubroutine compute_rec - 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. + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(kappa_rec_js), intent(in) :: self !< Kappa coefficient. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. -#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 + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu kappa coefficients object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) endfunction description elemental subroutine destroy(self) @@ -204,6 +204,6 @@ elemental subroutine destroy(self) class(kappa_rec_js), intent(inout) :: self !< Kappa. call self%destroy_ - if (allocated(self%values_rank_2)) deallocate(self%values_rank_2) + 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 index af182b6..d31df46 100644 --- a/src/lib/concrete_objects/wenoof_reconstructor_js.F90 +++ b/src/lib/concrete_objects/wenoof_reconstructor_js.F90 @@ -8,12 +8,12 @@ module wenoof_reconstructor_js #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_base_object -use wenoof_interpolations_factory -use wenoof_interpolations_object -use wenoof_interpolator_object -use wenoof_weights_factory -use wenoof_weights_object +use wenoof_base_object, only : base_object_constructor +use wenoof_interpolations_factory, only : interpolations_factory +use wenoof_interpolations_object, only : interpolations_object +use wenoof_interpolator_object, only : interpolator_object, interpolator_object_constructor +use wenoof_weights_factory, only : weights_factory +use wenoof_weights_object, only : weights_object implicit none private @@ -34,13 +34,13 @@ module wenoof_reconstructor_js !< 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_stencil_rank_1_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_stencil_rank_2_standard !< Interpolate values (without providing debug values). - procedure, pass(self) :: interpolate_stencil_rank_1_debug !< Interpolate values (providing also debug values). - procedure, pass(self) :: interpolate_stencil_rank_2_debug !< Interpolate values (providing also debug values). + procedure, pass(self) :: create !< Create reconstructor. + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy reconstructor. + procedure, pass(self) :: interpolate_int_debug !< Interpolate values (providing also debug values, interpolate). + procedure, pass(self) :: interpolate_int_standard !< Interpolate values (without providing debug values, interpolate). + procedure, pass(self) :: interpolate_rec_debug !< Interpolate values (providing also debug values, reconstruct). + procedure, pass(self) :: interpolate_rec_standard !< Interpolate values (without providing debug values, reconstruct). endtype reconstructor_js contains @@ -61,15 +61,18 @@ subroutine create(self, constructor) 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. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = 'Jiang-Shu reconstructor:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//self%weights%description() + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(reconstructor_js), intent(in) :: self !< Interpolator. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line character. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu reconstructor:'//NL + string = prefix_//string//' - S = '//trim(str(self%S))//NL + string = prefix_//string//self%weights%description(prefix=prefix_//' ') endfunction description elemental subroutine destroy(self) @@ -81,57 +84,62 @@ elemental subroutine destroy(self) if (allocated(self%weights)) deallocate(self%weights) endsubroutine destroy - pure subroutine interpolate_stencil_rank_1_debug(self, stencil, interpolation, si, weights) + pure subroutine interpolate_int_debug(self, stencil, interpolation, si, weights) + !< Interpolate values (providing also debug values, interpolate). + class(reconstructor_js), intent(in) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. + real(RPP), intent(out) :: interpolation !< Result of the interpolation. + real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [0:S-1]. + real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [0:S-1]. + ! empty procedure + endsubroutine interpolate_int_debug + + pure subroutine interpolate_int_standard(self, stencil, interpolation) + !< Interpolate values (without providing debug values, interpolate). + class(reconstructor_js), intent(in) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. + real(RPP), intent(out) :: interpolation !< Result of the interpolation. + ! empty procedure + endsubroutine interpolate_int_standard + + pure subroutine interpolate_rec_debug(self, stencil, interpolation, si, weights) !< Interpolate values (providing also debug values). - class(reconstructor_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. - real(RPP), intent(out) :: interpolation !< Result of the interpolation. - real(RPP), intent(out) :: si(0:) !< Computed values of smoothness indicators [0:S-1]. - real(RPP), intent(out) :: weights(0:) !< Weights of the stencils, [0:S-1]. - - ! Empty subroutine. - endsubroutine interpolate_stencil_rank_1_debug + !< @TODO implement smoothness indicator return. + class(reconstructor_js), intent(in) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. + real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. + real(RPP) :: interpolations(1:2, 0:self%S - 1) !< Stencils interpolations. + integer(I_P) :: f, s !< Counters. - pure subroutine interpolate_stencil_rank_2_debug(self, stencil, interpolation, si, weights) - !< Interpolate values (providing also debug values). - class(reconstructor_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(RPP), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1]. - real(RPP), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1]. - - call self%interpolate(stencil=stencil, interpolation=interpolation) - call self%weights%smoothness_indicators_of_rank_2(si=si) - !si = self%weights%smoothness_indicators() - ! weights = self%weights%values_rank_2 - endsubroutine interpolate_stencil_rank_2_debug - - pure subroutine interpolate_stencil_rank_1_standard(self, stencil, interpolation) - !< Interpolate values (without providing debug values). - class(reconstructor_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1 - self%S:) !< Stencil of the interpolation [1-S:-1+S]. - real(RPP), intent(out) :: interpolation !< Result of the interpolation. - - ! Empty subroutine. - endsubroutine interpolate_stencil_rank_1_standard + call self%interpolations%compute(stencil=stencil, values=interpolations) + call self%weights%compute(stencil=stencil, values=weights) + ! call self%weights%smoothness_indicators_of_rank_2(si=si) + interpolation = 0._RPP + do s=0, self%S - 1 ! stencils loop + do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) + interpolation(f) = interpolation(f) + weights(f, s) * interpolations(f, s) + enddo + enddo + endsubroutine interpolate_rec_debug - pure subroutine interpolate_stencil_rank_2_standard(self, stencil, interpolation) + pure subroutine interpolate_rec_standard(self, stencil, interpolation) !< Interpolate values (without providing debug values). - class(reconstructor_js), intent(inout) :: self !< Reconstructor. - real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. - real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. - real(RPP) :: interpolations(1:2, 0:self%S - 1) !< Stencils interpolations. - real(RPP) :: weights(1:2, 0:self%S - 1) !< Weights of stencils interpolations. - integer(I_P) :: f, s !< Counters. + class(reconstructor_js), intent(in) :: self !< Reconstructor. + real(RPP), intent(in) :: stencil(1:, 1 - self%S:) !< Stencil of the interpolation [1:2, 1-S:-1+S]. + real(RPP), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2]. + real(RPP) :: interpolations(1:2, 0:self%S - 1) !< Stencils interpolations. + real(RPP) :: weights(1:2, 0:self%S - 1) !< Weights of stencils interpolations. + integer(I_P) :: f, s !< Counters. call self%interpolations%compute(stencil=stencil, values=interpolations) call self%weights%compute(stencil=stencil, values=weights) interpolation = 0._RPP do s=0, self%S - 1 ! stencils loop do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - ! interpolation(f) = interpolation(f) + self%weights%values_rank_2(f, s) * self%interpolations%values_rank_2(f, s) interpolation(f) = interpolation(f) + weights(f, s) * interpolations(f, s) enddo enddo - endsubroutine interpolate_stencil_rank_2_standard + endsubroutine interpolate_rec_standard endmodule wenoof_reconstructor_js diff --git a/src/lib/concrete_objects/wenoof_weights_int_js.F90 b/src/lib/concrete_objects/wenoof_weights_int_js.F90 index 2de2047..bc0f586 100644 --- a/src/lib/concrete_objects/wenoof_weights_int_js.F90 +++ b/src/lib/concrete_objects/wenoof_weights_int_js.F90 @@ -12,19 +12,15 @@ module wenoof_weights_int_js #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_alpha_factory -use wenoof_alpha_object -use wenoof_alpha_int_js -use wenoof_alpha_int_m -use wenoof_alpha_int_z -use wenoof_base_object -use wenoof_beta_factory -use wenoof_beta_object -use wenoof_beta_int_js -use wenoof_kappa_factory -use wenoof_kappa_object -use wenoof_kappa_int_js -use wenoof_weights_object +use wenoof_alpha_factory, only : alpha_factory +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_base_object, only : base_object_constructor +use wenoof_beta_factory, only : beta_factory +use wenoof_beta_object, only : beta_object, beta_object_constructor +use wenoof_kappa_factory, only : kappa_factory +use wenoof_kappa_int_js, only : kappa_int_js +use wenoof_kappa_object, only : kappa_object, kappa_object_constructor +use wenoof_weights_object, only : weights_object, weights_object_constructor implicit none private @@ -50,13 +46,13 @@ module wenoof_weights_int_js class(kappa_object), allocatable :: kappa !< kappa coefficients (optimal, linear weights). contains ! deferred public methods - procedure, pass(self) :: create !< Create weights. - procedure, pass(self) :: compute_stencil_rank_1 !< Compute weights. - procedure, pass(self) :: compute_stencil_rank_2 !< Compute weights. - procedure, pass(self) :: description !< Return weights string-description. - procedure, pass(self) :: destroy !< Destroy weights. - procedure, pass(self) :: smoothness_indicators_of_rank_1 !< Return smoothness indicators. - procedure, pass(self) :: smoothness_indicators_of_rank_2 !< Return smoothness indicators. + procedure, pass(self) :: create !< Create weights. + procedure, pass(self) :: compute_int !< Compute weights (interpolate). + procedure, pass(self) :: compute_rec !< Compute weights (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy weights. + procedure, pass(self) :: smoothness_indicators_int !< Return smoothness indicators (interpolate). + procedure, pass(self) :: smoothness_indicators_rec !< Return smoothness indicators (reconstrcut). endtype weights_int_js contains @@ -71,72 +67,59 @@ subroutine create(self, constructor) call self%destroy call self%create_(constructor=constructor) - allocate(self%values_rank_1(0:self%S - 1)) - self%values_rank_1 = 0._RPP select type(constructor) type is(weights_int_js_constructor) associate(alpha_constructor=>constructor%alpha_constructor, & beta_constructor=>constructor%beta_constructor, & kappa_constructor=>constructor%kappa_constructor) - call a_factory%create(constructor=alpha_constructor, object=self%alpha) - ! select type(alpha_constructor) - ! type is(alpha_rec_js_constructor) - ! call factory%create(constructor=alpha_constructor, object=self%alpha) - ! type is(alpha_rec_m_constructor) - ! call factory%create(constructor=alpha_constructor, object=self%alpha) - ! type is(alpha_rec_z_constructor) - ! call factory%create(constructor=alpha_constructor, object=self%alpha) - ! endselect - call b_factory%create(constructor=beta_constructor, object=self%beta) - ! select type(beta_constructor) - ! type is(beta_rec_js_constructor) - ! allocate(beta_rec_js :: self%beta) - ! call self%beta%create(constructor=beta_constructor) - ! endselect - call k_factory%create(constructor=kappa_constructor, object=self%kappa) - ! 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_stencil_rank_1(self, stencil) + pure subroutine compute_int(self, stencil, values) !< Compute weights. - class(weights_int_js), intent(inout) :: self !< Weights. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. - integer(I_P) :: s !< Counters. - - call self%beta%compute(stencil=stencil) - call self%alpha%compute(beta=self%beta, kappa=self%kappa) + class(weights_int_js), intent(in) :: self !< Weights. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Weights values. + real(RPP) :: alpha(0:self%S-1) !< Aplha values. + real(RPP) :: beta(0:self%S-1) !< Beta values. + real(RPP) :: alpha_sum !< Sum of aplha values. + integer(I_P) :: s !< Counters. + + call self%beta%compute(stencil=stencil, values=beta) + select type(kappa => self%kappa) + class is(kappa_int_js) + call self%alpha%compute(beta=beta, kappa=kappa%values, values=alpha) + endselect + alpha_sum = sum(alpha) do s=0, self%S - 1 ! stencils loop - self%values_rank_1(s) = self%alpha%values_rank_1(s) / self%alpha%values_sum_rank_1 + values(s) = alpha(s) / alpha_sum enddo - endsubroutine compute_stencil_rank_1 + endsubroutine compute_int - pure subroutine compute_stencil_rank_2(self, stencil, values) + pure subroutine compute_rec(self, stencil, values) !< Compute weights. class(weights_int_js), intent(in) :: self !< Weights. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:,0:) !< Weights values of stencil interpolations. - - ! Empty routine. - endsubroutine compute_stencil_rank_2 - - pure function description(self) result(string) - !< Return string-description of weights. - class(weights_int_js), intent(in) :: self !< Weights. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Jiang-Shu weights:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//self%alpha%description() + ! empty procedure + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(weights_int_js), intent(in) :: self !< Weights. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu weights object for interpolation:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) + string = prefix_//string//self%alpha%description(prefix=prefix_//' ') endfunction description elemental subroutine destroy(self) @@ -144,29 +127,22 @@ elemental subroutine destroy(self) class(weights_int_js), intent(inout) :: self !< Weights. call self%destroy_ - if (allocated(self%values_rank_1)) deallocate(self%values_rank_1) if (allocated(self%alpha)) deallocate(self%alpha) if (allocated(self%beta)) deallocate(self%beta) if (allocated(self%kappa)) deallocate(self%kappa) endsubroutine destroy - pure subroutine smoothness_indicators_of_rank_1(self, si) - !< Return smoothness indicators.. + pure subroutine smoothness_indicators_int(self, si) + !< Return smoothness indicators (interpolate). class(weights_int_js), intent(in) :: self !< Weights. real(RPP), intent(out) :: si(:) !< Smoothness indicators. + ! TODO implement this + endsubroutine smoothness_indicators_int - if (allocated(self%beta)) then - if (allocated(self%beta%values_rank_1)) then - si = self%beta%values_rank_1 - endif - endif - endsubroutine smoothness_indicators_of_rank_1 - - pure subroutine smoothness_indicators_of_rank_2(self, si) - !< Return smoothness indicators.. + pure subroutine smoothness_indicators_rec(self, si) + !< Return smoothness indicators (reconstruct). class(weights_int_js), intent(in) :: self !< Weights. real(RPP), intent(out) :: si(:,:) !< Smoothness indicators. - - ! Empty routine - endsubroutine smoothness_indicators_of_rank_2 + ! empty procedure + endsubroutine smoothness_indicators_rec endmodule wenoof_weights_int_js diff --git a/src/lib/concrete_objects/wenoof_weights_rec_js.F90 b/src/lib/concrete_objects/wenoof_weights_rec_js.F90 index 7191b50..e8880b7 100644 --- a/src/lib/concrete_objects/wenoof_weights_rec_js.F90 +++ b/src/lib/concrete_objects/wenoof_weights_rec_js.F90 @@ -12,19 +12,15 @@ module wenoof_weights_rec_js #else use penf, only: I_P, RPP=>R8P, str #endif -use wenoof_alpha_factory -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_factory -use wenoof_beta_object -use wenoof_beta_rec_js -use wenoof_kappa_factory -use wenoof_kappa_object -use wenoof_kappa_rec_js -use wenoof_weights_object +use wenoof_alpha_factory, only : alpha_factory +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_base_object, only : base_object_constructor +use wenoof_beta_factory, only : beta_factory +use wenoof_beta_object, only : beta_object, beta_object_constructor +use wenoof_kappa_factory, only : kappa_factory +use wenoof_kappa_object, only : kappa_object, kappa_object_constructor +use wenoof_kappa_rec_js, only : kappa_rec_js +use wenoof_weights_object, only : weights_object, weights_object_constructor implicit none private @@ -50,13 +46,13 @@ module wenoof_weights_rec_js class(kappa_object), allocatable :: kappa !< kappa coefficients (optimal, linear weights). contains ! deferred public methods - procedure, pass(self) :: create !< Create weights. - procedure, pass(self) :: compute_stencil_rank_1 !< Compute weights. - procedure, pass(self) :: compute_stencil_rank_2 !< Compute weights. - procedure, pass(self) :: description !< Return weights string-description. - procedure, pass(self) :: destroy !< Destroy weights. - procedure, pass(self) :: smoothness_indicators_of_rank_1 !< Return smoothness indicators. - procedure, pass(self) :: smoothness_indicators_of_rank_2 !< Return smoothness indicators. + procedure, pass(self) :: create !< Create weights. + procedure, pass(self) :: compute_int !< Compute weights (interpolate). + procedure, pass(self) :: compute_rec !< Compute weights (reconstruct). + procedure, pass(self) :: description !< Return object string-description. + procedure, pass(self) :: destroy !< Destroy weights. + procedure, pass(self) :: smoothness_indicators_int !< Return smoothness indicators (interpolate). + procedure, pass(self) :: smoothness_indicators_rec !< Return smoothness indicators (reconstruct). endtype weights_rec_js contains @@ -76,68 +72,57 @@ subroutine create(self, constructor) associate(alpha_constructor=>constructor%alpha_constructor, & beta_constructor=>constructor%beta_constructor, & kappa_constructor=>constructor%kappa_constructor) - call a_factory%create(constructor=alpha_constructor, object=self%alpha) - ! select type(alpha_constructor) - ! type is(alpha_rec_js_constructor) - ! call factory%create(constructor=alpha_constructor, object=self%alpha) - ! type is(alpha_rec_m_constructor) - ! call factory%create(constructor=alpha_constructor, object=self%alpha) - ! type is(alpha_rec_z_constructor) - ! call factory%create(constructor=alpha_constructor, object=self%alpha) - ! endselect - call b_factory%create(constructor=beta_constructor, object=self%beta) - ! select type(beta_constructor) - ! type is(beta_rec_js_constructor) - ! allocate(beta_rec_js :: self%beta) - ! call self%beta%create(constructor=beta_constructor) - ! endselect - call k_factory%create(constructor=kappa_constructor, object=self%kappa) - ! 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_stencil_rank_1(self, stencil) - !< Compute weights. - class(weights_rec_js), intent(inout) :: self !< Weights. - real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + pure subroutine compute_int(self, stencil, values) + !< Compute weights (interpolate). + class(weights_rec_js), intent(in) :: self !< Weights. + real(RPP), intent(in) :: stencil(1-self%S:) !< Stencil used for the interpolation, [1-S:-1+S]. + real(RPP), intent(out) :: values(0:) !< Weights values. + ! empty procedure + endsubroutine compute_int - ! Empty routine. - endsubroutine compute_stencil_rank_1 - - pure subroutine compute_stencil_rank_2(self, stencil, values) - !< Compute weights. + pure subroutine compute_rec(self, stencil, values) + !< Compute weights (reconstruct). class(weights_rec_js), intent(in) :: self !< Weights. real(RPP), intent(in) :: stencil(1:,1-self%S:) !< Stencil used for the interpolation, [1:2, 1-S:-1+S]. real(RPP), intent(out) :: values(1:,0:) !< Weights values of stencil interpolations. + real(RPP) :: alpha(1:2,0:self%S-1) !< Aplha values. real(RPP) :: beta(1:2,0:self%S-1) !< Beta values. + real(RPP) :: alpha_sum(1:2) !< Sum of aplha values. integer(I_P) :: f, s !< Counters. call self%beta%compute(stencil=stencil, values=beta) - ! call self%alpha%compute(beta=self%beta, kappa=self%kappa) + select type(kappa => self%kappa) + class is(kappa_rec_js) + call self%alpha%compute(beta=beta, kappa=kappa%values, values=alpha) + endselect + alpha_sum(1) = sum(alpha(1,:)) + alpha_sum(2) = sum(alpha(2,:)) do s=0, self%S - 1 ! stencils loop do f=1, 2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) - values(f, s) = self%alpha%values_rank_2(f, s) / self%alpha%values_sum_rank_2(f) + values(f, s) = alpha(f, s) / alpha_sum(f) enddo enddo - endsubroutine compute_stencil_rank_2 - - pure function description(self) result(string) - !< Return string-description of weights. - class(weights_rec_js), intent(in) :: self !< Weights. - character(len=:), allocatable :: string !< String-description. - character(len=1), parameter :: nl=new_line('a') !< New line char. - - string = ' Jiang-Shu weights:'//nl - string = string//' - S = '//trim(str(self%S))//nl - string = string//self%alpha%description() + endsubroutine compute_rec + + pure function description(self, prefix) result(string) + !< Return object string-descripition. + class(weights_rec_js), intent(in) :: self !< Weights. + character(len=*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: string !< String-description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line char. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + string = prefix_//'Jiang-Shu weights object for reconstruction:'//NL + string = prefix_//string//' - S = '//trim(str(self%S)) + string = prefix_//string//self%alpha%description(prefix=prefix_//' ') endfunction description elemental subroutine destroy(self) @@ -150,23 +135,17 @@ elemental subroutine destroy(self) if (allocated(self%kappa)) deallocate(self%kappa) endsubroutine destroy - pure subroutine smoothness_indicators_of_rank_1(self, si) - !< Return smoothness indicators.. + pure subroutine smoothness_indicators_int(self, si) + !< Return smoothness indicators (interpolate). class(weights_rec_js), intent(in) :: self !< Weights. real(RPP), intent(out) :: si(:) !< Smoothness indicators. + ! empty procedure + endsubroutine smoothness_indicators_int - ! Empty routine - endsubroutine smoothness_indicators_of_rank_1 - - pure subroutine smoothness_indicators_of_rank_2(self, si) - !< Return smoothness indicators.. + pure subroutine smoothness_indicators_rec(self, si) + !< Return smoothness indicators (reconstruct). class(weights_rec_js), intent(in) :: self !< Weights. real(RPP), intent(out) :: si(:,:) !< Smoothness indicators. - - if (allocated(self%beta)) then - if (allocated(self%beta%values_rank_2)) then - si = self%beta%values_rank_2 - endif - endif - endsubroutine smoothness_indicators_of_rank_2 + ! TODO implement this + endsubroutine smoothness_indicators_rec endmodule wenoof_weights_rec_js diff --git a/src/lib/factories/wenoof_alpha_factory.f90 b/src/lib/factories/wenoof_alpha_factory.F90 similarity index 87% rename from src/lib/factories/wenoof_alpha_factory.f90 rename to src/lib/factories/wenoof_alpha_factory.F90 index 6835ef9..b05f0b7 100644 --- a/src/lib/factories/wenoof_alpha_factory.f90 +++ b/src/lib/factories/wenoof_alpha_factory.F90 @@ -8,13 +8,13 @@ module wenoof_alpha_factory #else use penf, only: I_P, RPP=>R8P #endif -use wenoof_alpha_object -use wenoof_alpha_rec_js -use wenoof_alpha_rec_m -use wenoof_alpha_rec_z -use wenoof_alpha_int_js -use wenoof_alpha_int_m -use wenoof_alpha_int_z +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_alpha_int_js, only : alpha_int_js, alpha_int_js_constructor +use wenoof_alpha_int_m, only : alpha_int_m, alpha_int_m_constructor +use wenoof_alpha_int_z, only : alpha_int_z, alpha_int_z_constructor +use wenoof_alpha_rec_js, only : alpha_rec_js, alpha_rec_js_constructor +use wenoof_alpha_rec_m, only : alpha_rec_m, alpha_rec_m_constructor +use wenoof_alpha_rec_z, only : alpha_rec_z, alpha_rec_z_constructor implicit none private @@ -93,7 +93,7 @@ subroutine create_constructor(interpolator_type, S, constructor, eps) endselect case('reconstructor-Z') allocate(alpha_rec_z_constructor :: constructor) - case default + case default write(stderr, '(A)') 'error: interpolator type "'//trim(adjustl(interpolator_type))//'" is unknown!' stop endselect diff --git a/src/lib/factories/wenoof_beta_factory.f90 b/src/lib/factories/wenoof_beta_factory.f90 index bbd399d..a93b505 100644 --- a/src/lib/factories/wenoof_beta_factory.f90 +++ b/src/lib/factories/wenoof_beta_factory.f90 @@ -2,10 +2,11 @@ module wenoof_beta_factory !< Wenoof beta factory. -use penf, only: I_P -use wenoof_beta_object -use wenoof_beta_rec_js -use wenoof_beta_int_js +use, intrinsic :: iso_fortran_env, only : stderr=>error_unit +use penf, only : I_P +use wenoof_beta_object, only : beta_object, beta_object_constructor +use wenoof_beta_rec_js, only : beta_rec_js, beta_rec_js_constructor +use wenoof_beta_int_js, only : beta_int_js, beta_int_js_constructor implicit none private @@ -59,6 +60,9 @@ subroutine create_constructor(interpolator_type, S, constructor) allocate(beta_rec_js_constructor :: constructor) case('reconstructor-Z') allocate(beta_rec_js_constructor :: constructor) + case default + write(stderr, '(A)') 'error: interpolator type "'//trim(adjustl(interpolator_type))//'" is unknown!' + stop endselect call constructor%create(S=S) endsubroutine create_constructor diff --git a/src/lib/factories/wenoof_interpolations_factory.f90 b/src/lib/factories/wenoof_interpolations_factory.F90 similarity index 91% rename from src/lib/factories/wenoof_interpolations_factory.f90 rename to src/lib/factories/wenoof_interpolations_factory.F90 index 1f36047..6487397 100644 --- a/src/lib/factories/wenoof_interpolations_factory.f90 +++ b/src/lib/factories/wenoof_interpolations_factory.F90 @@ -7,9 +7,9 @@ module wenoof_interpolations_factory #else use penf, only: I_P, RPP=>R8P #endif -use wenoof_interpolations_object -use wenoof_interpolations_rec_js -use wenoof_interpolations_int_js +use wenoof_interpolations_object, only : interpolations_object, interpolations_object_constructor +use wenoof_interpolations_rec_js, only : interpolations_rec_js, interpolations_rec_js_constructor +use wenoof_interpolations_int_js, only : interpolations_int_js, interpolations_int_js_constructor implicit none private diff --git a/src/lib/factories/wenoof_interpolator_factory.f90 b/src/lib/factories/wenoof_interpolator_factory.f90 index 7e66558..9606ffc 100644 --- a/src/lib/factories/wenoof_interpolator_factory.f90 +++ b/src/lib/factories/wenoof_interpolator_factory.f90 @@ -2,12 +2,13 @@ module wenoof_interpolator_factory !< Wenoof interpolator factory. +use, intrinsic :: iso_fortran_env, only : stderr=>error_unit use penf, only: I_P -use wenoof_interpolations_object -use wenoof_interpolator_object -use wenoof_interpolator_js -use wenoof_reconstructor_js -use wenoof_weights_object +use wenoof_interpolations_object, only : interpolations_object_constructor +use wenoof_interpolator_object, only : interpolator_object, interpolator_object_constructor +use wenoof_interpolator_js, only : interpolator_js, interpolator_js_constructor +use wenoof_reconstructor_js, only : reconstructor_js, reconstructor_js_constructor +use wenoof_weights_object, only : weights_object_constructor implicit none private @@ -41,29 +42,20 @@ subroutine create(constructor, object) subroutine create_constructor(interpolator_type, S, interpolations_constructor, weights_constructor, & constructor) !< Create an instance of concrete extension of [[weights_object_constructor]]. - character(*), intent(in) :: interpolator_type !< Type of interpolator. - integer(I_P), intent(in) :: S !< Stencils dimension. - class(interpolations_object_constructor), intent(in) :: interpolations_constructor !< Interpolations const. - class(weights_object_constructor), intent(in) :: weights_constructor !< Weights constructor. - class(interpolator_object_constructor), allocatable, intent(out) :: constructor !< Constructor. + character(*), intent(in) :: interpolator_type !< Type of interpolator. + integer(I_P), intent(in) :: S !< Stencils dimension. + class(interpolations_object_constructor), intent(in) :: interpolations_constructor !< Interpolations const. + class(weights_object_constructor), intent(in) :: weights_constructor !< Weights constructor. + class(interpolator_object_constructor), allocatable, intent(out) :: constructor !< Constructor. select case(trim(adjustl(interpolator_type))) - case('interpolator-JS') + case('interpolator-JS', 'interpolator-M-JS', 'interpolator-M-Z', 'interpolator-Z') allocate(interpolator_js_constructor :: constructor) - case('interpolator-M-JS') - allocate(interpolator_js_constructor :: constructor) - case('interpolator-M-Z') - allocate(interpolator_js_constructor :: constructor) - case('interpolator-Z') - allocate(interpolator_js_constructor :: constructor) - case('reconstructor-JS') - allocate(reconstructor_js_constructor :: constructor) - case('reconstructor-M-JS') - allocate(reconstructor_js_constructor :: constructor) - case('reconstructor-M-Z') - allocate(reconstructor_js_constructor :: constructor) - case('reconstructor-Z') + case('reconstructor-JS', 'reconstructor-M-JS', 'reconstructor-M-Z', 'reconstructor-Z') allocate(reconstructor_js_constructor :: constructor) + case default + write(stderr, '(A)') 'error: interpolator type "'//trim(adjustl(interpolator_type))//'" is unknown!' + stop endselect call constructor%create(S=S) select type(constructor) diff --git a/src/lib/factories/wenoof_kappa_factory.f90 b/src/lib/factories/wenoof_kappa_factory.f90 index 27e803b..c155485 100644 --- a/src/lib/factories/wenoof_kappa_factory.f90 +++ b/src/lib/factories/wenoof_kappa_factory.f90 @@ -7,10 +7,10 @@ module wenoof_kappa_factory #else use penf, only: I_P, RPP=>R8P #endif -use wenoof_interpolations_object -use wenoof_kappa_object -use wenoof_kappa_rec_js -use wenoof_kappa_int_js +use wenoof_interpolations_object, only : interpolations_object_constructor +use wenoof_kappa_object, only : kappa_object, kappa_object_constructor +use wenoof_kappa_rec_js, only : kappa_rec_js, kappa_rec_js_constructor +use wenoof_kappa_int_js, only : kappa_int_js, kappa_int_js_constructor implicit none private @@ -21,10 +21,10 @@ module wenoof_kappa_factory contains ! public methods procedure, nopass :: create !< Create a concrete instance of [[kappa_object]]. - procedure, nopass :: create_constructor_rec - procedure, nopass :: create_constructor_int - generic :: create_constructor => create_constructor_rec, & !< Create a concrete instance - create_constructor_int !< of [[kappa_object_constructor]]. + procedure, nopass :: create_constructor_int !< Create [[kappa_object_constructor]] (interpolate). + procedure, nopass :: create_constructor_rec !< Create [[kappa_object_constructor]] (reconstruct). + generic :: create_constructor => create_constructor_int, & !< Create a concrete instance + create_constructor_rec !< of [[kappa_object_constructor]]. endtype kappa_factory contains @@ -64,12 +64,11 @@ subroutine create_constructor_int(interpolator_type, S, stencil, x_target, inter class(kappa_object_constructor), allocatable, intent(out) :: constructor !< Constructor. allocate(kappa_int_js_constructor :: constructor) - allocate(constructor%stencil(1-S:S-1)) - constructor%stencil = stencil - constructor%x_target = x_target - call constructor%create(S=S) select type(constructor) type is(kappa_int_js_constructor) + constructor%stencil = stencil + constructor%x_target = x_target + call constructor%create(S=S) allocate(constructor%interpolations_constructor, source=interpolations_constructor) endselect endsubroutine create_constructor_int diff --git a/src/lib/factories/wenoof_objects_factory.f90 b/src/lib/factories/wenoof_objects_factory.f90 index 58007c7..1a18873 100644 --- a/src/lib/factories/wenoof_objects_factory.f90 +++ b/src/lib/factories/wenoof_objects_factory.f90 @@ -7,18 +7,18 @@ module wenoof_objects_factory #else use penf, only: I_P, RPP=>R8P #endif -use wenoof_alpha_factory -use wenoof_alpha_object -use wenoof_beta_factory -use wenoof_beta_object -use wenoof_kappa_factory -use wenoof_kappa_object -use wenoof_interpolations_factory -use wenoof_interpolations_object -use wenoof_interpolator_factory -use wenoof_interpolator_object -use wenoof_weights_factory -use wenoof_weights_object +use wenoof_alpha_factory, only : alpha_factory +use wenoof_alpha_object, only : alpha_object, alpha_object_constructor +use wenoof_beta_factory, only : beta_factory +use wenoof_beta_object, only : beta_object, beta_object_constructor +use wenoof_kappa_factory, only : kappa_factory +use wenoof_kappa_object, only : kappa_object, kappa_object_constructor +use wenoof_interpolations_factory, only : interpolations_factory +use wenoof_interpolations_object, only : interpolations_object, interpolations_object_constructor +use wenoof_interpolator_factory, only : interpolator_factory +use wenoof_interpolator_object, only : interpolator_object, interpolator_object_constructor +use wenoof_weights_factory, only : weights_factory +use wenoof_weights_object, only : weights_object, weights_object_constructor implicit none private diff --git a/src/lib/factories/wenoof_weights_factory.f90 b/src/lib/factories/wenoof_weights_factory.f90 index c339060..89df666 100644 --- a/src/lib/factories/wenoof_weights_factory.f90 +++ b/src/lib/factories/wenoof_weights_factory.f90 @@ -2,13 +2,14 @@ module wenoof_weights_factory !< Wenoof weights factory. +use, intrinsic :: iso_fortran_env, only : stderr=>error_unit use penf, only: I_P -use wenoof_alpha_object -use wenoof_beta_object -use wenoof_kappa_object -use wenoof_weights_object -use wenoof_weights_int_js -use wenoof_weights_rec_js +use wenoof_alpha_object , only : alpha_object_constructor +use wenoof_beta_object , only : beta_object_constructor +use wenoof_kappa_object , only : kappa_object_constructor +use wenoof_weights_object, only : weights_object, weights_object_constructor +use wenoof_weights_int_js, only : weights_int_js, weights_int_js_constructor +use wenoof_weights_rec_js, only : weights_rec_js, weights_rec_js_constructor implicit none private @@ -66,6 +67,9 @@ subroutine create_constructor(interpolator_type, S, alpha_constructor, beta_cons allocate(weights_rec_js_constructor :: constructor) case('reconstructor-Z') allocate(weights_rec_js_constructor :: constructor) + case default + write(stderr, '(A)') 'error: interpolator type "'//trim(adjustl(interpolator_type))//'" is unknown!' + stop endselect call constructor%create(S=S) select type(constructor)