Changes to support H rescaling

src/ALE/P1M_functions.F90

@@ -41,7 +41,7 @@ module P1M_functions
 !------------------------------------------------------------------------------
 ! p1m interpolation
 !------------------------------------------------------------------------------
-subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients )
+subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect )
 ! ------------------------------------------------------------------------------
 ! Linearly interpolate between edge values.
 ! The resulting piecewise interpolant is stored in 'ppoly'.
@@ -57,18 +57,23 @@ subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients )
 ! ------------------------------------------------------------------------------
 
   ! Arguments
-  integer,            intent(in)    :: N ! Number of cells
-  real, dimension(:), intent(in)    :: h ! cell widths (size N)
-  real, dimension(:), intent(in)    :: u ! cell averages (size N)
-  real, dimension(:,:), intent(inout) :: ppoly_E
-  real, dimension(:,:), intent(inout) :: ppoly_coefficients
+  integer,              intent(in)    :: N !< Number of cells
+  real, dimension(:),   intent(in)    :: h !< cell widths (size N)
+  real, dimension(:),   intent(in)    :: u !< cell average properties (size N)
+  real, dimension(:,:), intent(inout) :: ppoly_E !< Potentially modified edge values,
+                                           !! with the same units as u.
+  real, dimension(:,:), intent(inout) :: ppoly_coefficients !< Potentially modified
+                                           !! piecewise polynomial coefficients, mainly
+                                           !! with the same units as u.
+  real,       optional, intent(in)    :: h_neglect !< A negligibly small width
+                                           !! in the same units as h.
 
   ! Local variables
   integer   :: k            ! loop index
   real      :: u0_l, u0_r   ! edge values (left and right)
 
   ! Bound edge values (routine found in 'edge_values.F90')
-  call bound_edge_values( N, h, u, ppoly_E )
+  call bound_edge_values( N, h, u, ppoly_E, h_neglect )
 
   ! Systematically average discontinuous edge values (routine found in
   ! 'edge_values.F90')

src/ALE/P3M_functions.F90

@@ -22,14 +22,16 @@ module P3M_functions
 public P3M_interpolation
 public P3M_boundary_extrapolation
 
-real, parameter :: h_neglect = 1.E-30
+real, parameter :: hNeglect_dflt = 1.E-30
+real, parameter :: hNeglect_edge_dflt = 1.E-10
 
 contains
 
 !------------------------------------------------------------------------------
 ! p3m interpolation
 ! -----------------------------------------------------------------------------
-subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
+subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, &
+                              h_neglect )
 !------------------------------------------------------------------------------
 ! Cubic interpolation between edges.
 !
@@ -47,23 +49,25 @@ subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
   real, dimension(:,:), intent(inout) :: ppoly_E            !Edge value of polynomial
   real, dimension(:,:), intent(inout) :: ppoly_S            !Edge slope of polynomial
   real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial
-
+  real,       optional, intent(in)    :: h_neglect !< A negligibly small width for the
+                                          !! purpose of cell reconstructions
+                                          !! in the same units as h.
 
   ! Call the limiter for p3m, which takes care of everything from
   ! computing the coefficients of the cubic to monotonizing it.
   ! This routine could be called directly instead of having to call
   ! 'P3M_interpolation' first but we do that to provide an homogeneous
   ! interface.
 
-  call P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
+  call P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect )
 
 end subroutine P3M_interpolation
 
 
 !------------------------------------------------------------------------------
 ! p3m limiter
 ! -----------------------------------------------------------------------------
-subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
+subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect )
 !------------------------------------------------------------------------------
 ! The p3m limiter operates as follows:
 !
@@ -84,25 +88,29 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
   real, dimension(:,:), intent(inout) :: ppoly_E            !Edge value of polynomial
   real, dimension(:,:), intent(inout) :: ppoly_S            !Edge slope of polynomial
   real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial
-
-!  real, dimension(:,:), intent(inout) :: ppoly_coefficients
+  real,       optional, intent(in)    :: h_neglect !< A negligibly small width for
+                                           !! the purpose of cell reconstructions
+                                           !! in the same units as h.
 
   ! Local variables
-  integer   :: k            ! loop index
-  integer   :: monotonic    ! boolean indicating whether the cubic is monotonic
-  real      :: u0_l, u0_r   ! edge values
-  real      :: u1_l, u1_r   ! edge slopes
-  real      :: u_l, u_c, u_r        ! left, center and right cell averages
-  real      :: h_l, h_c, h_r        ! left, center and right cell widths
-  real      :: sigma_l, sigma_c, sigma_r    ! left, center and right
-                                            ! van Leer slopes
-  real      :: slope        ! retained PLM slope
-  real      :: eps
+  integer :: k            ! loop index
+  integer :: monotonic    ! boolean indicating whether the cubic is monotonic
+  real    :: u0_l, u0_r   ! edge values
+  real    :: u1_l, u1_r   ! edge slopes
+  real    :: u_l, u_c, u_r        ! left, center and right cell averages
+  real    :: h_l, h_c, h_r        ! left, center and right cell widths
+  real    :: sigma_l, sigma_c, sigma_r    ! left, center and right
+                                          ! van Leer slopes
+  real    :: slope        ! retained PLM slope
+  real    :: eps
+  real    :: hNeglect
+
+  hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect
 
   eps = 1e-10
 
   ! 1. Bound edge values (boundary cells are assumed to be local extrema)
-  call bound_edge_values( N, h, u, ppoly_E )
+  call bound_edge_values( N, h, u, ppoly_E, hNeglect )
 
   ! 2. Systematically average discontinuous edge values
   call average_discontinuous_edge_values( N, ppoly_E )
@@ -142,9 +150,9 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
     end if
 
     ! Compute limited slope
-    sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + h_neglect )
-    sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + h_neglect )
-    sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + h_neglect )
+    sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hNeglect )
+    sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hNeglect )
+    sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hNeglect )
 
     if ( (sigma_l * sigma_r) .GT. 0.0 ) then
       slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
@@ -154,12 +162,12 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
 
     ! If the slopes are close to zero in machine precision and in absolute
     ! value, we set the slope to zero. This prevents asymmetric representation
-    ! near extrema.
-    if ( abs(u1_l*h_c) .LT. eps ) then
+    ! near extrema.  These expressions are both nondimensional.
+    if ( abs(u1_l*h_c) < eps ) then
       u1_l = 0.0
     end if
 
-    if ( abs(u1_r*h_c) .LT. eps ) then
+    if ( abs(u1_r*h_c) < eps ) then
       u1_r = 0.0
     end if
 
@@ -201,7 +209,8 @@ end subroutine P3M_limiter
 !------------------------------------------------------------------------------
 ! p3m boundary extrapolation
 ! -----------------------------------------------------------------------------
-subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients )
+subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, &
+                                       h_neglect, h_neglect_edge )
 !------------------------------------------------------------------------------
 ! The following explanations apply to the left boundary cell. The same
 ! reasoning holds for the right boundary cell.
@@ -222,25 +231,33 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici
   real, dimension(:,:), intent(inout) :: ppoly_E            !Edge value of polynomial
   real, dimension(:,:), intent(inout) :: ppoly_S            !Edge slope of polynomial
   real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial
+  real,       optional, intent(in)    :: h_neglect !< A negligibly small width for the
+                                          !! purpose of cell reconstructions
+                                          !! in the same units as h.
+  real,       optional, intent(in)    :: h_neglect_edge !< A negligibly small width
+                                          !! for the purpose of finding edge values
+                                          !! in the same units as h.
 
   ! Local variables
-  integer       :: i0, i1
-  integer       :: monotonic
-  real          :: u0, u1
-  real          :: h0, h1
-  real          :: b, c, d
-  real          :: u0_l, u0_r
-  real          :: u1_l, u1_r
-  real          :: eps
-  real          :: slope
-
-  eps = 1e-10
+  integer :: i0, i1
+  integer :: monotonic
+  real    :: u0, u1
+  real    :: h0, h1
+  real    :: b, c, d
+  real    :: u0_l, u0_r
+  real    :: u1_l, u1_r
+  real    :: slope
+  real    :: hNeglect, hNeglect_edge
+
+  hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect
+  hNeglect_edge = hNeglect_edge_dflt
+  if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge
 
   ! ----- Left boundary -----
   i0 = 1
   i1 = 2
-  h0 = h(i0) + eps
-  h1 = h(i1) + eps
+  h0 = h(i0) + hNeglect_edge
+  h1 = h(i1) + hNeglect_edge
   u0 = u(i0)
   u1 = u(i1)
 
@@ -250,7 +267,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici
   u1_r = b / h1     ! derivative evaluated at xi = 0.0, expressed w.r.t. x
 
   ! Limit the right slope by the PLM limited slope
-  slope = 2.0 * ( u1 - u0 ) / ( h0 + h_neglect )
+  slope = 2.0 * ( u1 - u0 ) / ( h0 + hNeglect )
   if ( abs(u1_r) .GT. abs(slope) ) then
     u1_r = slope
   end if
@@ -263,7 +280,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici
   ! edge value and slope by computing the parabola as determined by
   ! the right edge value and slope and the boundary cell average
   u0_l = 3.0 * u0 + 0.5 * h0*u1_r - 2.0 * u0_r
-  u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + h_neglect )
+  u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + hNeglect )
 
   ! Check whether the edge values are monotonic. For example, if the left edge
   ! value is larger than the right edge value while the slope is positive, the
@@ -297,8 +314,8 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici
   ! ----- Right boundary -----
   i0 = N-1
   i1 = N
-  h0 = h(i0) + eps
-  h1 = h(i1) + eps
+  h0 = h(i0) + hNeglect_edge
+  h1 = h(i1) + hNeglect_edge
   u0 = u(i0)
   u1 = u(i1)
 
@@ -307,10 +324,10 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici
   b = ppoly_coefficients(i0,2)
   c = ppoly_coefficients(i0,3)
   d = ppoly_coefficients(i0,4)
-  u1_l = (b + 2*c + 3*d) / ( h0 + h_neglect ) ! derivative evaluated at xi = 1.0
+  u1_l = (b + 2*c + 3*d) / ( h0 + hNeglect ) ! derivative evaluated at xi = 1.0
 
   ! Limit the left slope by the PLM limited slope
-  slope = 2.0 * ( u1 - u0 ) / ( h1 + h_neglect )
+  slope = 2.0 * ( u1 - u0 ) / ( h1 + hNeglect )
   if ( abs(u1_l) .GT. abs(slope) ) then
     u1_l = slope
   end if
@@ -323,7 +340,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici
   ! edge value and slope by computing the parabola as determined by
   ! the left edge value and slope and the boundary cell average
   u0_r = 3.0 * u1 - 0.5 * h1*u1_l - 2.0 * u0_l
-  u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + h_neglect )
+  u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + hNeglect )
 
   ! Check whether the edge values are monotonic. For example, if the right edge
   ! value is smaller than the left edge value while the slope is positive, the

src/ALE/PLM_functions.F90

@@ -16,14 +16,14 @@ module PLM_functions
 
 public PLM_reconstruction, PLM_boundary_extrapolation
 
-real, parameter :: h_neglect = 1.E-30
+real, parameter :: hNeglect_dflt = 1.E-30
 
 contains
 
 !------------------------------------------------------------------------------
 ! PLM_reconstruction
 ! -----------------------------------------------------------------------------
-subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients )
+subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect )
 !------------------------------------------------------------------------------
 ! Reconstruction by linear polynomials within each cell.
 !
@@ -43,6 +43,9 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients )
   real, dimension(:),   intent(in)    :: u ! cell averages (size N)
   real, dimension(:,:), intent(inout) :: ppoly_E
   real, dimension(:,:), intent(inout) :: ppoly_coefficients
+  real,       optional, intent(in)    :: h_neglect !< A negligibly small width for
+                                          !! the purpose of cell reconstructions
+                                          !! in the same units as h
 
   ! Local variables
   integer       :: k                    ! loop index
@@ -55,6 +58,9 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients )
   real          :: u_min, u_max, e_l, e_r, edge
   real          :: almost_one, almost_two
   real, dimension(N) :: slp, mslp
+  real    :: hNeglect
+
+  hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect
 
   almost_one = 1. - epsilon(slope)
   almost_two = 2. * almost_one
@@ -67,7 +73,7 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients )
 
     ! Get cell widths
     h_l = h(k-1) ; h_c = h(k) ; h_r = h(k+1)
-    h_cn = max( h_c, h_neglect ) ! To avoid division by zero
+    h_cn = max( h_c, hNeglect ) ! To avoid division by zero
 
     ! Side differences
     sigma_r = u_r - u_c
@@ -83,7 +89,7 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients )
 
     ! This is the original estimate of the second order slope from Laurent
     ! but multiplied by h_c
-    sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + h_neglect) )
+    sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + hNeglect) )
 
     if ( (sigma_l * sigma_r) > 0.0 ) then
       ! This limits the slope so that the edge values are bounded by the
@@ -209,7 +215,7 @@ end subroutine PLM_reconstruction
 !------------------------------------------------------------------------------
 ! plm boundary extrapolation
 ! -----------------------------------------------------------------------------
-subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients )
+subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect )
 !------------------------------------------------------------------------------
 ! Reconstruction by linear polynomials within boundary cells.
 ! The left and right edge values in the left and right boundary cells,
@@ -233,17 +239,23 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients )
   real, dimension(:),   intent(in)    :: u ! cell averages (size N)
   real, dimension(:,:), intent(inout) :: ppoly_E
   real, dimension(:,:), intent(inout) :: ppoly_coefficients
+  real,       optional, intent(in)    :: h_neglect !< A negligibly small width for
+                                          !! the purpose of cell reconstructions
+                                          !! in the same units as h
 
   ! Local variables
-  real          :: u0, u1               ! cell averages
-  real          :: h0, h1               ! corresponding cell widths
-  real          :: slope                ! retained PLM slope
+  real    :: u0, u1               ! cell averages
+  real    :: h0, h1               ! corresponding cell widths
+  real    :: slope                ! retained PLM slope
+  real    :: hNeglect
+
+  hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect
 
   ! -----------------------------------------
   ! Left edge value in the left boundary cell
   ! -----------------------------------------
-  h0 = h(1) + h_neglect
-  h1 = h(2) + h_neglect
+  h0 = h(1) + hNeglect
+  h1 = h(2) + hNeglect
 
   u0 = u(1)
   u1 = u(2)
@@ -264,8 +276,8 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients )
   ! ------------------------------------------
   ! Right edge value in the left boundary cell
   ! ------------------------------------------
-  h0 = h(N-1) + h_neglect
-  h1 = h(N) + h_neglect
+  h0 = h(N-1) + hNeglect
+  h1 = h(N) + hNeglect
 
   u0 = u(N-1)
   u1 = u(N)