Skip to content

Commit

Permalink
add weights export
Browse files Browse the repository at this point in the history
  • Loading branch information
@szaghi
szaghi committed Jan 16, 2017
1 parent 66a370f commit 9175e22
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 42 deletions.
5 changes: 3 additions & 2 deletions src/lib/wenoof_interpolator.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,14 +104,15 @@ pure subroutine interpolate_standard(self, S, stencil, location, interpolation)
#endif
endsubroutine interpolate_standard

pure subroutine interpolate_debug(self, S, stencil, location, interpolation, si)
!< Interpolate values (without providing debug values).
pure subroutine interpolate_debug(self, S, stencil, location, interpolation, si, weights)
!< Interpolate values (providing also debug values).
class(interpolator), intent(inout) :: self !< Interpolator.
integer(I_P), intent(in) :: S !< Number of stencils actually used.
real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil of the interpolation [1:2, 1-S:-1+S].
character(*), intent(in) :: location !< Location of interpolation: left, right, both.
real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2].
real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1].
real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1].

#ifndef DEBUG
! error stop in pure procedure is a F2015 feature not yet supported in debug mode
Expand Down
105 changes: 71 additions & 34 deletions src/lib/wenoof_interpolator_js.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@ module wenoof_interpolator_js
! public methods
procedure, pass(self) :: create !< Create interpolator.
procedure, pass(self) :: destroy !< Destroy interpolator.
! private methods
procedure, pass(self), private :: compute_convolution !< Compute convolution.
procedure, pass(self), private :: compute_weights !< Compute weights.
endtype interpolator_js

contains
Expand Down Expand Up @@ -116,58 +119,43 @@ pure subroutine interpolate_standard(self, S, stencil, location, interpolation)
integer(I_P) :: f1, f2, ff !< Faces to be computed.
integer(I_P) :: f, k !< Counters.

select case(location)
case('both', 'b')
f1=1_I_P; f2=2_I_P; ff=0_I_P
case('left', 'l')
f1=1_I_P; f2=1_I_P; ff=0_I_P
case('right', 'r')
f1=2_I_P; f2=2_I_P; ff=-1_I_P
endselect
call compute_faces_indexes(location=location, f1=f1, f2=f2, ff=ff)

call self%polynom%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff)

call self%is%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff)

call self%alpha%compute(S=S, weight_opt=self%weights%opt, IS = self%IS%si, eps = self%eps, f1=f1, f2=f2)
call self%alpha%compute(S=S, weight_opt=self%weights%opt, IS=self%IS%si, eps=self%eps, f1=f1, f2=f2)

! computing the weights
do k = 0, S - 1 ! stencils loop
do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2)
weights(f, k) = self%alpha%alpha_coef(f, k) / self%alpha%alpha_tot(f)
enddo
enddo
call self%compute_weights(S=S, f1=f1, f2=f2, ff=ff, weights=weights)

! computing the convolution
interpolation = 0._R_P
do k = 0, S - 1 ! stencils loop
do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2)
interpolation(f + ff) = interpolation(f + ff) + weights(f, k) * self%polynom%poly(f, k)
enddo
enddo
call self%compute_convolution(S=S, f1=f1, f2=f2, ff=ff, weights=weights, interpolation=interpolation)
endsubroutine interpolate_standard

pure subroutine interpolate_debug(self, S, stencil, location, interpolation, si)
!< Interpolate values (without providing debug values).
pure subroutine interpolate_debug(self, S, stencil, location, interpolation, si, weights)
!< Interpolate values (providing also debug values).
class(interpolator_js), intent(inout) :: self !< Interpolator.
integer(I_P), intent(in) :: S !< Number of stencils actually used.
real(R_P), intent(in) :: stencil(1:, 1 - S:) !< Stencil of the interpolation [1:2, 1-S:-1+S].
character(*), intent(in) :: location !< Location of interpolation: left, right, both.
real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2].
real(R_P), intent(out) :: si(1:, 0:) !< Computed values of smoothness indicators [1:2, 0:S-1].
real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1].
integer(I_P) :: f1, f2, ff !< Faces to be computed.
integer(I_P) :: f !< Counter.
integer(I_P) :: f, k !< Counters.

select case(location)
case('both', 'b')
f1=1_I_P; f2=2_I_P; ff=0_I_P
case('left', 'l')
f1=1_I_P; f2=1_I_P; ff=0_I_P
case('right', 'r')
f1=2_I_P; f2=2_I_P; ff=-1_I_P
endselect
call compute_faces_indexes(location=location, f1=f1, f2=f2, ff=ff)

call self%polynom%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff)

call self%is%compute(S=S, stencil=stencil, f1=f1, f2=f2, ff = ff)

call self%alpha%compute(S=S, weight_opt=self%weights%opt, IS=self%IS%si, eps=self%eps, f1=f1, f2=f2)

call self%compute_weights(S=S, f1=f1, f2=f2, ff=ff, weights=weights)

call self%compute_convolution(S=S, f1=f1, f2=f2, ff=ff, weights=weights, interpolation=interpolation)

call self%interpolate_standard(S=S, stencil=stencil, location=location, interpolation=interpolation)
associate(is => self%is)
select type(is)
class is(smoothness_indicators)
Expand Down Expand Up @@ -201,4 +189,53 @@ elemental subroutine destroy(self)
self%S = 0_I_P
self%eps = 0._R_P
endsubroutine destroy

! private methods
pure subroutine compute_convolution(self, S, f1, f2, ff, weights, interpolation)
!< Compute convolution.
class(interpolator_js), intent(in) :: self !< Interpolator.
integer(I_P), intent(in) :: S !< Number of stencils actually used.
integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed.
real(R_P), intent(in) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1].
real(R_P), intent(out) :: interpolation(1:) !< Result of the interpolation, [1:2].
integer(I_P) :: f, k !< Counters.

interpolation = 0._R_P
do k = 0, S - 1 ! stencils loop
do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2)
interpolation(f + ff) = interpolation(f + ff) + weights(f + ff, k) * self%polynom%poly(f, k)
enddo
enddo
endsubroutine compute_convolution

pure subroutine compute_weights(self, S, f1, f2, ff, weights)
!< Compute weights.
class(interpolator_js), intent(in) :: self !< Interpolator.
integer(I_P), intent(in) :: S !< Number of stencils actually used.
integer(I_P), intent(in) :: f1, f2, ff !< Faces to be computed.
real(R_P), intent(out) :: weights(1:, 0:) !< Weights of the stencils, [1:2, 0:S-1].
integer(I_P) :: f, k !< Counters.

do k = 0, S - 1 ! stencils loop
do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2)
weights(f + ff, k) = self%alpha%alpha_coef(f, k) / self%alpha%alpha_tot(f)
enddo
enddo
endsubroutine compute_weights

! private non TBP
pure subroutine compute_faces_indexes(location, f1, f2, ff)
!< Compute faces indexes given the queried location.
character(*), intent(in) :: location !< Location of interpolation: left, right, both.
integer(I_P), intent(out) :: f1, f2, ff !< Faces to be computed.

select case(location)
case('both', 'b')
f1=1_I_P; f2=2_I_P; ff=0_I_P
case('left', 'l')
f1=1_I_P; f2=1_I_P; ff=0_I_P
case('right', 'r')
f1=2_I_P; f2=2_I_P; ff=-1_I_P
endselect
endsubroutine compute_faces_indexes
endmodule wenoof_interpolator_js
20 changes: 14 additions & 6 deletions src/tests/sin_reconstruction.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ module test_module
real(R_P), allocatable :: fx_face(:) !< Face reference values [1:points_number].
real(R_P), allocatable :: interpolation(:) !< Interpolated values [1:points_number].
real(R_P), allocatable :: si(:,:) !< Computed smoothness indicators [1:points_number,0:S-1].
real(R_P), allocatable :: weights(:,:) !< Computed weights [1:points_number,0:S-1].
real(R_P) :: error_L2 !< L2 norm of the numerical error.
endtype solution_data

Expand Down Expand Up @@ -105,12 +106,14 @@ subroutine allocate_solution_data(self)
allocate(self%solution(pn, s)%fx_face( 1:self%points_number(pn) ))
allocate(self%solution(pn, s)%interpolation( 1:self%points_number(pn) ))
allocate(self%solution(pn, s)%si( 1:self%points_number(pn), 0:self%S(s)-1))
allocate(self%solution(pn, s)%weights( 1:self%points_number(pn), 0:self%S(s)-1))
self%solution(pn, s)%x_cell = 0._R_P
self%solution(pn, s)%fx_cell = 0._R_P
self%solution(pn, s)%x_face = 0._R_P
self%solution(pn, s)%fx_face = 0._R_P
self%solution(pn, s)%interpolation = 0._R_P
self%solution(pn, s)%si = 0._R_P
self%solution(pn, s)%weights = 0._R_P
enddo
enddo
endsubroutine allocate_solution_data
Expand Down Expand Up @@ -249,7 +252,8 @@ subroutine perform(self)
shape=[1,2*self%S(s)-1]), &
location='right', &
interpolation=self%solution(pn, s)%interpolation(i:i), &
si=self%solution(pn, s)%si(i:i, 0:self%S(s)-1))
si=self%solution(pn, s)%si(i:i, 0:self%S(s)-1), &
weights=self%solution(pn, s)%weights(i:i, 0:self%S(s)-1))
enddo
enddo
enddo
Expand Down Expand Up @@ -283,15 +287,19 @@ subroutine save_results_and_plots(self)
do ss=0, self%S(s)-1
buffer = buffer//' "si-'//trim(str(ss, .true.))//'"'
enddo
do ss=0, self%S(s)-1
buffer = buffer//' "W-'//trim(str(ss, .true.))//'"'
enddo
write(file_unit, "(A)") buffer
write(file_unit, "(A)") 'ZONE T = "'//'S_'//trim(str(self%S(s), .true.))//&
'-Np_'//trim(str(self%points_number(pn), .true.))//'"'
do i = 1, self%points_number(pn)
write(file_unit, "("//trim(str(3+self%S(s), .true.))//"("//FR_P//",1X))") &
self%solution(pn, s)%x_face(i), &
self%solution(pn, s)%fx_face(i), &
self%solution(pn, s)%interpolation(i), &
(self%solution(pn, s)%si(i, ss), ss=0, self%S(s)-1)
write(file_unit, "("//trim(str(3+2*self%S(s), .true.))//"("//FR_P//",1X))") &
self%solution(pn, s)%x_face(i), &
self%solution(pn, s)%fx_face(i), &
self%solution(pn, s)%interpolation(i), &
(self%solution(pn, s)%si(i, ss), ss=0, self%S(s)-1), &
(self%solution(pn, s)%weights(i, ss), ss=0, self%S(s)-1)
enddo
close(file_unit)
enddo
Expand Down

0 comments on commit 9175e22

Please sign in to comment.