Fortran-FOSS-Programmers / WenOOF

WENO interpolation Object Oriented Fortran library
35 stars 17 forks source link

Add methods to save computed IS for debugging aims #17

Closed szaghi closed 7 years ago

szaghi commented 7 years ago

For debugging aims it can be very useful to save the computed values of smoothness indicators. Add it to the interpolate methods, better with dynamic dispatch overloading, e.g. interpolate method becomes a generic

type, extends(base_object) :: interpolator
  !< Abstract interpolator object.
  !<
  !< @note Do not implement any actual interpolator: provide the interface for the different interpolators implemented.
  class(smoothness_indicators), allocatable :: is      !< Smoothness indicators.
  class(alpha_coefficients),    allocatable :: alpha   !< Alpha coefficients.
  class(optimal_weights),       allocatable :: weights !< Optimal weights.
  class(polynomials),           allocatable :: polynom !< Polynomilas.
  contains
    ! public deferred methods
    procedure, nopass     :: description !< Return interpolator string-description.
    procedure, pass(self) :: interpolate_standard !< Interpolate values, without providing debug values.    
    procedure, pass(self) :: interpolate_debug !< Interpolate values, providing also debug values.
    ! public methods
    generic :: interpolate => interpolate_standard, interpolate_debug !< Interpolate values.
    procedure, pass(self) :: create  !< Create interpolator.
    procedure, pass(self) :: destroy !< Destroy interpolator.
endtype interpolator

contains

  pure subroutine interpolate_standard(self, S, stencil, location, interpolation)
  !< Interpolate values.
  class(interpolator), intent(inout) :: self                  !< Interpolator.
  integer(I_P),        intent(in)    :: S                     !< Number of stencils actually used.
  real(R_P),           intent(in)    :: stencil(1:, 1 - S:)   !< Stencil of the interpolation [1:2, 1-S:-1+S].
  character(*),        intent(in)    :: location              !< Location of interpolation: left, right, both.
  real(R_P),           intent(out)   :: interpolation(1:)     !< Result of the interpolation, [1:2].

#ifndef DEBUG
  ! error stop in pure procedure is a F2015 feature not yet supported in debug mode
  error stop 'interpolator%interpolate to be implemented by your concrete interpolator object'
#endif
  endsubroutine interpolate_standard

  pure subroutine interpolate_debug(self, S, stencil, location, interpolation, smoothness_indicators, ...)
  !< Interpolate 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(in)    :: smoothness_indicatorsl(1:, 1 - S:)   !< Stencil of the interpolation [1:2, 1-S:-1+S].
...

#ifndef DEBUG
  ! error stop in pure procedure is a F2015 feature not yet supported in debug mode
  error stop 'interpolator%interpolate to be implemented by your concrete interpolator object'
#endif
  endsubroutine interpolate_debug
zbeekman commented 7 years ago

It can also be useful if you're trying to come up with non-linear optimizations (i.e., adaptation limiters etc.) Martin has a paper on this (Martin & Wu I think... should be somewhere here: http://croccolab.umd.edu/publications/journals.html) where the adaptation (smoothness indicators) are visualized for shock wave/BL interactions.

szaghi commented 7 years ago

@zbeekman

It can also be useful if you're trying to come up with non-linear optimizations

Sure, non-linear optimizations will be investigated.

where the adaptation (smoothness indicators) are visualized for shock wave/BL interactions.

Yep, I want a debug mode just for research/debug on schemes and their possible modification. Wave/boundary-layer is one of my interest (I suspect they play a role in SRM): do you know the work of Pirozzoli?

Cheers

szaghi commented 7 years ago

done.