From 5e33043c10b799973c9700fa23a3839a8067167f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 6 May 2018 09:30:35 -0400 Subject: [PATCH] Added dOxyGen comments in ALE regridding code Added dOxyGen comments for all of the subroutines and their arguments in the PXM_functions and regrid_... modules in src/ALE. Also shortened the name of several variables. All answers are bitwise identical. --- src/ALE/P1M_functions.F90 | 38 +++---- src/ALE/P3M_functions.F90 | 149 +++++++++++++----------- src/ALE/PCM_functions.F90 | 21 ++-- src/ALE/PLM_functions.F90 | 70 ++++++------ src/ALE/PPM_functions.F90 | 57 +++++----- src/ALE/PQM_functions.F90 | 146 ++++++++++++------------ src/ALE/polynomial_functions.F90 | 87 +++++++------- src/ALE/regrid_edge_slopes.F90 | 32 +++--- src/ALE/regrid_edge_values.F90 | 190 ++++++++++++++----------------- src/ALE/regrid_interp.F90 | 24 ++-- src/ALE/regrid_solvers.F90 | 27 ++--- 11 files changed, 420 insertions(+), 421 deletions(-) diff --git a/src/ALE/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index aafaec2580..8590a7297f 100644 --- a/src/ALE/P1M_functions.F90 +++ b/src/ALE/P1M_functions.F90 @@ -39,9 +39,8 @@ module P1M_functions !------------------------------------------------------------------------------ -! p1m interpolation -!------------------------------------------------------------------------------ -subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) +!> Linearly interpolate between edge values +subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coef, h_neglect ) ! ------------------------------------------------------------------------------ ! Linearly interpolate between edge values. ! The resulting piecewise interpolant is stored in 'ppoly'. @@ -62,7 +61,7 @@ subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) 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 + real, dimension(:,:), intent(inout) :: ppoly_coef !< Potentially modified !! piecewise polynomial coefficients, mainly !! with the same units as u. real, optional, intent(in) :: h_neglect !< A negligibly small width @@ -85,8 +84,8 @@ subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) u0_l = ppoly_E(k,1) u0_r = ppoly_E(k,2) - ppoly_coefficients(k,1) = u0_l - ppoly_coefficients(k,2) = u0_r - u0_l + ppoly_coef(k,1) = u0_l + ppoly_coef(k,2) = u0_r - u0_l end do ! end loop on interior cells @@ -94,9 +93,8 @@ end subroutine P1M_interpolation !------------------------------------------------------------------------------ -! p1m boundary extrapolation -! ----------------------------------------------------------------------------- -subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) +!> Interpolation by linear polynomials within boundary cells +subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) !------------------------------------------------------------------------------ ! Interpolation by linear polynomials within boundary cells. ! The left and right edge values in the left and right boundary cells, @@ -106,18 +104,20 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) ! h: thicknesses of grid cells ! u: cell averages to use in constructing piecewise polynomials ! ppoly_E : edge values of piecewise polynomials -! ppoly_coefficients : coefficients of piecewise polynomials +! ppoly_coef : coefficients of piecewise polynomials ! ! It is assumed that the size of the array 'u' is equal to the number of cells ! defining 'grid' and 'ppoly'. No consistency check is performed here. !------------------------------------------------------------------------------ ! 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 averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !< edge values of piecewise polynomials, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly + !! with the same units as u. ! Local variables real :: u0, u1 ! cell averages @@ -157,8 +157,8 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) ppoly_E(1,1) = u0 end if - ppoly_coefficients(1,1) = ppoly_E(1,1) - ppoly_coefficients(1,2) = ppoly_E(1,2) - ppoly_E(1,1) + ppoly_coef(1,1) = ppoly_E(1,1) + ppoly_coef(1,2) = ppoly_E(1,2) - ppoly_E(1,1) ! ------------------------------------------ ! Right edge value in the left boundary cell @@ -183,8 +183,8 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) ppoly_E(N,2) = u1 end if - ppoly_coefficients(N,1) = ppoly_E(N,1) - ppoly_coefficients(N,2) = ppoly_E(N,2) - ppoly_E(N,1) + ppoly_coef(N,1) = ppoly_E(N,1) + ppoly_coef(N,2) = ppoly_E(N,2) - ppoly_E(N,1) end subroutine P1M_boundary_extrapolation diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index a532ca7003..acc3e064ce 100644 --- a/src/ALE/P3M_functions.F90 +++ b/src/ALE/P3M_functions.F90 @@ -28,9 +28,9 @@ module P3M_functions contains !------------------------------------------------------------------------------ -! p3m interpolation -! ----------------------------------------------------------------------------- -subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, & +!> Set up a piecewise cubic cubic interpolation from cell averages and estimated +!! edge slopes and values +subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & h_neglect ) !------------------------------------------------------------------------------ ! Cubic interpolation between edges. @@ -43,12 +43,15 @@ subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, 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 !Edge value of polynomial - real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial - real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + 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 !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, + !! in the units of u over the units of h. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly + !! with the same units as u. real, optional, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions !! in the same units as h. @@ -59,15 +62,15 @@ subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, & ! 'P3M_interpolation' first but we do that to provide an homogeneous ! interface. - call P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect ) + call P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) end subroutine P3M_interpolation !------------------------------------------------------------------------------ -! p3m limiter -! ----------------------------------------------------------------------------- -subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect ) +!> Adust a piecewise cubic reconstruction with a limiter that adjusts the edge +!! values and slopes +subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) !------------------------------------------------------------------------------ ! The p3m limiter operates as follows: ! @@ -82,12 +85,14 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect !------------------------------------------------------------------------------ ! 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 !Edge value of polynomial - real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial - real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + 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 !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, + !! in the units of u over the units of h. + real, dimension(:,:), intent(inout) :: ppoly_coef !< 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. @@ -182,10 +187,10 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect end if ! Build cubic interpolant (compute the coefficients) - call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) + call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef ) ! Check whether cubic is monotonic - monotonic = is_cubic_monotonic( ppoly_coefficients, k ) + monotonic = is_cubic_monotonic( ppoly_coef, k ) ! If cubic is not monotonic, monotonize it by modifiying the ! edge slopes, store the new edge slopes and recompute the @@ -199,7 +204,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect ppoly_S(k,2) = u1_r ! Recompute coefficients of cubic - call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) + call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef ) end do ! loop on cells @@ -207,9 +212,9 @@ end subroutine P3M_limiter !------------------------------------------------------------------------------ -! p3m boundary extrapolation -! ----------------------------------------------------------------------------- -subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, & +!> calculate the edge values and slopes at boundary cells as part of building a +!! piecewise peicewise cubic sub-grid scale profiles +subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & h_neglect, h_neglect_edge ) !------------------------------------------------------------------------------ ! The following explanations apply to the left boundary cell. The same @@ -225,12 +230,15 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici !------------------------------------------------------------------------------ ! 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 !Edge value of polynomial - real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial - real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + 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 !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, + !! in the units of u over the units of h. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly + !! with the same units as u. real, optional, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions !! in the same units as h. @@ -263,7 +271,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly_coefficients(i1,2) + b = ppoly_coef(i1,2) u1_r = b / h1 ! derivative evaluated at xi = 0.0, expressed w.r.t. x ! Limit the right slope by the PLM limited slope @@ -298,8 +306,8 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici ppoly_S(i0,2) = u1_r ! Store edge values and slope, build cubic and check monotonicity - call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coefficients ) - monotonic = is_cubic_monotonic( ppoly_coefficients, i0 ) + call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coef ) + monotonic = is_cubic_monotonic( ppoly_coef, i0 ) if ( monotonic == 0 ) then call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r ) @@ -307,7 +315,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici ! Rebuild cubic after monotonization ppoly_S(i0,1) = u1_l ppoly_S(i0,2) = u1_r - call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coefficients ) + call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coef ) end if @@ -321,9 +329,9 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici ! Compute the right edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly_coefficients(i0,2) - c = ppoly_coefficients(i0,3) - d = ppoly_coefficients(i0,4) + b = ppoly_coef(i0,2) + c = ppoly_coef(i0,3) + d = ppoly_coef(i0,4) u1_l = (b + 2*c + 3*d) / ( h0 + hNeglect ) ! derivative evaluated at xi = 1.0 ! Limit the left slope by the PLM limited slope @@ -357,8 +365,8 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici ppoly_S(i1,1) = u1_l ppoly_S(i1,2) = u1_r - call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coefficients ) - monotonic = is_cubic_monotonic( ppoly_coefficients, i1 ) + call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coef ) + monotonic = is_cubic_monotonic( ppoly_coef, i1 ) if ( monotonic == 0 ) then call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r ) @@ -366,7 +374,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici ! Rebuild cubic after monotonization ppoly_S(i1,1) = u1_l ppoly_S(i1,2) = u1_r - call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coefficients ) + call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coef ) end if @@ -374,9 +382,8 @@ end subroutine P3M_boundary_extrapolation !------------------------------------------------------------------------------ -! Build cubic interpolant in cell k -! ----------------------------------------------------------------------------- -subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) +!> Build cubic interpolant in cell k +subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef ) !------------------------------------------------------------------------------ ! Given edge values and edge slopes, compute coefficients of cubic in cell k. ! @@ -385,11 +392,14 @@ subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Arguments - real, dimension(:), intent(in) :: h ! cell widths (size N) - integer, intent(in) :: k - real, dimension(:,:), intent(in) :: ppoly_E !Edge value of polynomial - real, dimension(:,:), intent(in) :: ppoly_S !Edge slope of polynomial - real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + real, dimension(:), intent(in) :: h !< cell widths (size N) + integer, intent(in) :: k !< The index of the cell to work on + real, dimension(:,:), intent(in) :: ppoly_E !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(in) :: ppoly_S !< Edge slope of polynomial, + !! in the units of u over the units of h. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly + !! with the same units as u. ! Local variables real :: u0_l, u0_r ! edge values @@ -410,18 +420,17 @@ subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) a2 = 3.0 * ( u0_r - u0_l ) - u1_r - 2.0 * u1_l a3 = u1_r + u1_l + 2.0 * ( u0_l - u0_r ) - ppoly_coefficients(k,1) = a0 - ppoly_coefficients(k,2) = a1 - ppoly_coefficients(k,3) = a2 - ppoly_coefficients(k,4) = a3 + ppoly_coef(k,1) = a0 + ppoly_coef(k,2) = a1 + ppoly_coef(k,3) = a2 + ppoly_coef(k,4) = a3 end subroutine build_cubic_interpolant !------------------------------------------------------------------------------ -! Check whether cubic is monotonic -! ----------------------------------------------------------------------------- -integer function is_cubic_monotonic( ppoly_coefficients, k ) +!> Check whether the cubic reconstruction in cell k is monotonic +integer function is_cubic_monotonic( ppoly_coef, k ) !------------------------------------------------------------------------------ ! This function checks whether the cubic curve in cell k is monotonic. ! If so, returns 1. Otherwise, returns 0. @@ -432,8 +441,8 @@ integer function is_cubic_monotonic( ppoly_coefficients, k ) !------------------------------------------------------------------------------ ! Arguments - real, dimension(:,:), intent(in) :: ppoly_coefficients - integer, intent(in) :: k + real, dimension(:,:), intent(in) :: ppoly_coef !< Coefficients of cubic polynomial + integer, intent(in) :: k !< The index of the cell to work on ! Local variables integer :: monotonic ! boolean indicating if monotonic or not @@ -447,10 +456,10 @@ integer function is_cubic_monotonic( ppoly_coefficients, k ) ! to be equal to 0 or 1, respectively eps = 1e-14 - a0 = ppoly_coefficients(k,1) - a1 = ppoly_coefficients(k,2) - a2 = ppoly_coefficients(k,3) - a3 = ppoly_coefficients(k,4) + a0 = ppoly_coef(k,1) + a1 = ppoly_coef(k,2) + a2 = ppoly_coef(k,3) + a3 = ppoly_coef(k,4) a = a1 b = 2.0 * a2 @@ -490,8 +499,7 @@ end function is_cubic_monotonic !------------------------------------------------------------------------------ -! Monotonize cubic curve -! ----------------------------------------------------------------------------- +!> Monotonize a cubic curve by modifying the edge slopes. subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ) !------------------------------------------------------------------------------ ! This routine takes care of monotonizing a cubic on [0,1] by modifying the @@ -522,11 +530,14 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r !------------------------------------------------------------------------------ ! Arguments - real, intent(in) :: h ! cell width - real, intent(in) :: u0_l, u0_r ! edge values - real, intent(in) :: sigma_l, sigma_r ! left and right 2nd-order slopes - real, intent(in) :: slope ! limited PLM slope - real, intent(inout) :: u1_l, u1_r ! edge slopes + real, intent(in) :: h !< cell width + real, intent(in) :: u0_l !< left edge value + real, intent(in) :: u0_r !< right edge value + real, intent(in) :: sigma_l !< left 2nd-order slopes + real, intent(in) :: sigma_r !< right 2nd-order slopes + real, intent(in) :: slope !< limited PLM slope + real, intent(inout) :: u1_l !< left edge slopes + real, intent(inout) :: u1_r !< right edge slopes ! Local variables integer :: found_ip diff --git a/src/ALE/PCM_functions.F90 b/src/ALE/PCM_functions.F90 index b09f6e080e..bcb963faa6 100644 --- a/src/ALE/PCM_functions.F90 +++ b/src/ALE/PCM_functions.F90 @@ -19,9 +19,10 @@ module PCM_functions contains !------------------------------------------------------------------------------ -! pcm_reconstruction -!------------------------------------------------------------------------------ -subroutine PCM_reconstruction( N, u, ppoly_E, ppoly_coefficients ) +!> Reconstruction by constant polynomials within each cell. There is nothing to +!! do but this routine is provided to ensure a homogeneous interface +!! throughout the regridding toolbox. +subroutine PCM_reconstruction( N, u, ppoly_E, ppoly_coef ) !------------------------------------------------------------------------------ ! Reconstruction by constant polynomials within each cell. There is nothing to ! do but this routine is provided to ensure a homogeneous interface @@ -31,24 +32,26 @@ subroutine PCM_reconstruction( N, u, ppoly_E, ppoly_coefficients ) ! h: thicknesses of grid cells ! u: cell averages to use in constructing piecewise polynomials ! ppoly_E : edge values of piecewise polynomials -! ppoly_coefficients : coefficients of piecewise polynomials +! ppoly_coef : coefficients of piecewise polynomials ! ! It is assumed that the dimension of 'u' is equal to the number of cells ! defining 'grid' and 'ppoly'. No consistency check is performed. !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: u ! cell averages - real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial - real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + integer, intent(in) :: N !< Number of cells + real, dimension(:), intent(in) :: u !< cell averages + real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, + !! with the same units as u. ! Local variables integer :: k ! The coefficients of the piecewise constant polynomial are simply ! the cell averages. - ppoly_coefficients(:,1) = u(:) + ppoly_coef(:,1) = u(:) ! The edge values are equal to the cell average do k = 1,N diff --git a/src/ALE/PLM_functions.F90 b/src/ALE/PLM_functions.F90 index 83eea1518b..73f9206c21 100644 --- a/src/ALE/PLM_functions.F90 +++ b/src/ALE/PLM_functions.F90 @@ -21,9 +21,8 @@ module PLM_functions contains !------------------------------------------------------------------------------ -! PLM_reconstruction -! ----------------------------------------------------------------------------- -subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) +!> Reconstruction by linear polynomials within each cell +subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coef, h_neglect ) !------------------------------------------------------------------------------ ! Reconstruction by linear polynomials within each cell. ! @@ -31,21 +30,23 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) ! h: thicknesses of grid cells ! u: cell averages to use in constructing piecewise polynomials ! ppoly_E : edge values of piecewise polynomials -! ppoly_coefficients : coefficients of piecewise polynomials +! ppoly_coef : coefficients of piecewise polynomials ! ! It is assumed that the size of the array 'u' is equal to the number of cells ! defining 'grid' and 'ppoly'. No consistency check is performed here. !------------------------------------------------------------------------------ ! 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 averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !< edge values of piecewise polynomials, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly + !! with the same units as u. real, optional, intent(in) :: h_neglect !< A negligibly small width for - !! the purpose of cell reconstructions - !! in the same units as h + !! the purpose of cell reconstructions + !! in the same units as h ! Local variables integer :: k ! loop index @@ -171,8 +172,8 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) ! Store and return edge values and polynomial coefficients. ppoly_E(1,1) = u(1) ppoly_E(1,2) = u(1) - ppoly_coefficients(1,1) = u(1) - ppoly_coefficients(1,2) = 0. + ppoly_coef(1,1) = u(1) + ppoly_coef(1,2) = 0. do k = 2, N-1 slope = sign( mslp(k), slp(k) ) u_l = u(k) - 0.5 * slope ! Left edge value of cell k @@ -194,28 +195,27 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) ppoly_E(k,1) = u_l ppoly_E(k,2) = u_r - ppoly_coefficients(k,1) = u_l - ppoly_coefficients(k,2) = ( u_r - u_l ) + ppoly_coef(k,1) = u_l + ppoly_coef(k,2) = ( u_r - u_l ) ! Check to see if this evaluation of the polynomial at x=1 would be ! monotonic w.r.t. the next cell's edge value. If not, scale back! - edge = ppoly_coefficients(k,2) + ppoly_coefficients(k,1) + edge = ppoly_coef(k,2) + ppoly_coef(k,1) e_r = u(k+1) - 0.5 * sign( mslp(k+1), slp(k+1) ) if ( (edge-u(k))*(e_r-edge)<0.) then - ppoly_coefficients(k,2) = ppoly_coefficients(k,2) * almost_one + ppoly_coef(k,2) = ppoly_coef(k,2) * almost_one endif enddo ppoly_E(N,1) = u(N) ppoly_E(N,2) = u(N) - ppoly_coefficients(N,1) = u(N) - ppoly_coefficients(N,2) = 0. + ppoly_coef(N,1) = u(N) + ppoly_coef(N,2) = 0. end subroutine PLM_reconstruction !------------------------------------------------------------------------------ -! plm boundary extrapolation -! ----------------------------------------------------------------------------- -subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect ) +!> Reconstruction by linear polynomials within boundary cells +subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef, h_neglect ) !------------------------------------------------------------------------------ ! Reconstruction by linear polynomials within boundary cells. ! The left and right edge values in the left and right boundary cells, @@ -227,21 +227,23 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n ! h: thicknesses of grid cells ! u: cell averages to use in constructing piecewise polynomials ! ppoly_E : edge values of piecewise polynomials -! ppoly_coefficients : coefficients of piecewise polynomials +! ppoly_coef : coefficients of piecewise polynomials ! ! It is assumed that the size of the array 'u' is equal to the number of cells ! defining 'grid' and 'ppoly'. No consistency check is performed here. !------------------------------------------------------------------------------ ! 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 averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !< edge values of piecewise polynomials, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly + !! with the same units as u. real, optional, intent(in) :: h_neglect !< A negligibly small width for - !! the purpose of cell reconstructions - !! in the same units as h + !! the purpose of cell reconstructions + !! in the same units as h ! Local variables real :: u0, u1 ! cell averages @@ -270,8 +272,8 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n ppoly_E(1,1) = u0 - 0.5 * slope ppoly_E(1,2) = u0 + 0.5 * slope - ppoly_coefficients(1,1) = ppoly_E(1,1) - ppoly_coefficients(1,2) = ppoly_E(1,2) - ppoly_E(1,1) + ppoly_coef(1,1) = ppoly_E(1,1) + ppoly_coef(1,2) = ppoly_E(1,2) - ppoly_E(1,1) ! ------------------------------------------ ! Right edge value in the left boundary cell @@ -292,8 +294,8 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n ppoly_E(N,1) = u1 - 0.5 * slope ppoly_E(N,2) = u1 + 0.5 * slope - ppoly_coefficients(N,1) = ppoly_E(N,1) - ppoly_coefficients(N,2) = ppoly_E(N,2) - ppoly_E(N,1) + ppoly_coef(N,1) = ppoly_E(N,1) + ppoly_coef(N,2) = ppoly_E(N,2) - ppoly_E(N,1) end subroutine PLM_boundary_extrapolation diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index 484e2e0d40..d0eb8325ad 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -25,12 +25,14 @@ module PPM_functions contains !> Builds quadratic polynomials coefficients from cell mean and edge values. -subroutine PPM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect) +subroutine PPM_reconstruction( N, h, u, ppoly_E, ppoly_coef, h_neglect) integer, intent(in) :: N !< Number of cells real, dimension(N), intent(in) :: h !< Cell widths real, dimension(N), intent(in) :: u !< Cell averages - real, dimension(N,2), intent(inout) :: ppoly_E !< Edge values - real, dimension(N,3), intent(inout) :: ppoly_coefficients !< Polynomial coefficients + real, dimension(N,2), intent(inout) :: ppoly_E !< Edge values, + !! with the same units as u. + real, dimension(N,3), intent(inout) :: ppoly_coef !< 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 @@ -47,9 +49,9 @@ subroutine PPM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients, h_neglect) edge_r = ppoly_E(k,2) ! Store polynomial coefficients - ppoly_coefficients(k,1) = edge_l - ppoly_coefficients(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r ) - ppoly_coefficients(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) ) + ppoly_coef(k,1) = edge_l + ppoly_coef(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r ) + ppoly_coef(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) ) enddo @@ -127,9 +129,8 @@ end subroutine PPM_limiter_standard !------------------------------------------------------------------------------ -! ppm boundary extrapolation -! ----------------------------------------------------------------------------- -subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_neglect) +!> Reconstruction by parabolas within boundary cells +subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef, h_neglect) !------------------------------------------------------------------------------ ! Reconstruction by parabolas within boundary cells. ! @@ -148,21 +149,23 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n ! h: thicknesses of grid cells ! u: cell averages to use in constructing piecewise polynomials ! ppoly_E : edge values of piecewise polynomials -! ppoly_coefficients : coefficients of piecewise polynomials +! ppoly_coef : coefficients of piecewise polynomials ! ! It is assumed that the size of the array 'u' is equal to the number of cells ! defining 'grid' and 'ppoly'. No consistency check is performed here. !------------------------------------------------------------------------------ ! 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 - real, optional, intent(in) :: h_neglect !< A negligibly small width for - !! the purpose of cell reconstructions - !! in the same units as h. + 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 !< edge values of piecewise polynomials, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly + !! with the same units as u. + 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 :: i0, i1 @@ -187,7 +190,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly_coefficients(i1,2) + b = ppoly_coef(i1,2) u1_r = b *((h0+hNeglect)/(h1+hNeglect)) ! derivative evaluated at xi = 0.0, ! expressed w.r.t. xi (local coord. system) @@ -225,9 +228,9 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u0 ) - ppoly_coefficients(i0,1) = a - ppoly_coefficients(i0,2) = b - ppoly_coefficients(i0,3) = c + ppoly_coef(i0,1) = a + ppoly_coef(i0,2) = b + ppoly_coef(i0,3) = c ! ----- Right boundary ----- i0 = N-1 @@ -239,8 +242,8 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n ! Compute the right edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly_coefficients(i0,2) - c = ppoly_coefficients(i0,3) + b = ppoly_coef(i0,2) + c = ppoly_coef(i0,3) u1_l = (b + 2*c) ! derivative evaluated at xi = 1.0 u1_l = u1_l * ((h1+hNeglect)/(h0+hNeglect)) @@ -278,9 +281,9 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients, h_n b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u1 ) - ppoly_coefficients(i1,1) = a - ppoly_coefficients(i1,2) = b - ppoly_coefficients(i1,3) = c + ppoly_coef(i1,1) = a + ppoly_coef(i1,2) = b + ppoly_coef(i1,3) = c end subroutine PPM_boundary_extrapolation diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index 8d15e6dd98..6c89c7ac10 100644 --- a/src/ALE/PQM_functions.F90 +++ b/src/ALE/PQM_functions.F90 @@ -22,9 +22,8 @@ module PQM_functions contains !------------------------------------------------------------------------------ -! PQM_reconstruction -! ----------------------------------------------------------------------------- -subroutine PQM_reconstruction( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect ) +!> PQM_reconstruction does reconstruction by quartic polynomials within each cell. +subroutine PQM_reconstruction( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) !------------------------------------------------------------------------------ ! Reconstruction by quartic polynomials within each cell. ! @@ -37,15 +36,18 @@ subroutine PQM_reconstruction( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_ !------------------------------------------------------------------------------ ! 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 !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 + 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 !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, + !! in the units of u over the units of h. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly + !! with the same units as u. + 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 @@ -75,11 +77,11 @@ subroutine PQM_reconstruction( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_ e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r) ! Store coefficients - ppoly_coefficients(k,1) = a - ppoly_coefficients(k,2) = b - ppoly_coefficients(k,3) = c - ppoly_coefficients(k,4) = d - ppoly_coefficients(k,5) = e + ppoly_coef(k,1) = a + ppoly_coef(k,2) = b + ppoly_coef(k,3) = c + ppoly_coef(k,4) = d + ppoly_coef(k,5) = e end do ! end loop on cells @@ -87,8 +89,7 @@ end subroutine PQM_reconstruction !------------------------------------------------------------------------------ -! Limit pqm -! ----------------------------------------------------------------------------- +!> Limit the piecewise quartic method reconstruction subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) !------------------------------------------------------------------------------ ! Standard PQM limiter (White & Adcroft, JCP 2008). @@ -369,9 +370,8 @@ end subroutine PQM_limiter !------------------------------------------------------------------------------ -! pqm boundary extrapolation -! ----------------------------------------------------------------------------- -subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) +!> piecewise quartic method boundary extrapolation +subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) !------------------------------------------------------------------------------ ! Reconstruction by parabolas within boundary cells. ! @@ -395,11 +395,13 @@ subroutine PQM_boundary_extrapolation( 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 !Edge value of polynomial - real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + 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 !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly + !! with the same units as u. ! Local variables integer :: i0, i1 @@ -421,7 +423,7 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly_coefficients(i1,2) + b = ppoly_coef(i1,2) u1_r = b *(h0/h1) ! derivative evaluated at xi = 0.0, ! expressed w.r.t. xi (local coord. system) @@ -460,11 +462,11 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) c = 3.0 * ( u0_r + u0_l - 2.0 * u0 ) ! The quartic is reduced to a parabola in the boundary cell - ppoly_coefficients(i0,1) = a - ppoly_coefficients(i0,2) = b - ppoly_coefficients(i0,3) = c - ppoly_coefficients(i0,4) = 0.0 - ppoly_coefficients(i0,5) = 0.0 + ppoly_coef(i0,1) = a + ppoly_coef(i0,2) = b + ppoly_coef(i0,3) = c + ppoly_coef(i0,4) = 0.0 + ppoly_coef(i0,5) = 0.0 ! ----- Right boundary ----- i0 = N-1 @@ -476,10 +478,10 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) ! Compute the right edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly_coefficients(i0,2) - c = ppoly_coefficients(i0,3) - d = ppoly_coefficients(i0,4) - e = ppoly_coefficients(i0,5) + b = ppoly_coef(i0,2) + c = ppoly_coef(i0,3) + d = ppoly_coef(i0,4) + e = ppoly_coef(i0,5) u1_l = (b + 2*c + 3*d + 4*e) ! derivative evaluated at xi = 1.0 u1_l = u1_l * (h1/h0) @@ -518,19 +520,18 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) c = 3.0 * ( u0_r + u0_l - 2.0 * u1 ) ! The quartic is reduced to a parabola in the boundary cell - ppoly_coefficients(i1,1) = a - ppoly_coefficients(i1,2) = b - ppoly_coefficients(i1,3) = c - ppoly_coefficients(i1,4) = 0.0 - ppoly_coefficients(i1,5) = 0.0 + ppoly_coef(i1,1) = a + ppoly_coef(i1,2) = b + ppoly_coef(i1,3) = c + ppoly_coef(i1,4) = 0.0 + ppoly_coef(i1,5) = 0.0 end subroutine PQM_boundary_extrapolation !------------------------------------------------------------------------------ -! pqm boundary extrapolation using rational function -! ----------------------------------------------------------------------------- -subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect ) +!> pqm boundary extrapolation using a rational function +subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) !------------------------------------------------------------------------------ ! Reconstruction by parabolas within boundary cells. ! @@ -554,15 +555,18 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff !------------------------------------------------------------------------------ ! 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 !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. + 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 !< Edge value of polynomial, + !! with the same units as u. + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, + !! in the units of u over the units of h. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly + !! with the same units as u. + 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 :: i0, i1 @@ -600,8 +604,8 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ! The right edge value and slope of the boundary cell are taken to be the ! left edge value and slope of the adjacent cell - a = ppoly_coefficients(i1,1) - b = ppoly_coefficients(i1,2) + a = ppoly_coef(i1,1) + b = ppoly_coef(i1,2) u0_r = a ! edge value u1_r = b / (h1 + hNeglect) ! edge slope (w.r.t. global coord.) @@ -725,11 +729,11 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r) ! Store coefficients - ppoly_coefficients(i0,1) = a - ppoly_coefficients(i0,2) = b - ppoly_coefficients(i0,3) = c - ppoly_coefficients(i0,4) = d - ppoly_coefficients(i0,5) = e + ppoly_coef(i0,1) = a + ppoly_coef(i0,2) = b + ppoly_coef(i0,3) = c + ppoly_coef(i0,4) = d + ppoly_coef(i0,5) = e ! ----- Right boundary (BOTTOM) ----- i0 = N-1 @@ -747,11 +751,11 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ! The left edge value and slope of the boundary cell are taken to be the ! right edge value and slope of the adjacent cell - a = ppoly_coefficients(i0,1) - b = ppoly_coefficients(i0,2) - c = ppoly_coefficients(i0,3) - d = ppoly_coefficients(i0,4) - e = ppoly_coefficients(i0,5) + a = ppoly_coef(i0,1) + b = ppoly_coef(i0,2) + c = ppoly_coef(i0,3) + d = ppoly_coef(i0,4) + e = ppoly_coef(i0,5) u0_l = a + b + c + d + e ! edge value u1_l = (b + 2*c + 3*d + 4*e) / h0 ! edge slope (w.r.t. global coord.) @@ -877,11 +881,11 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff d = -60.0 * um + h1 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r) - ppoly_coefficients(i1,1) = a - ppoly_coefficients(i1,2) = b - ppoly_coefficients(i1,3) = c - ppoly_coefficients(i1,4) = d - ppoly_coefficients(i1,5) = e + ppoly_coef(i1,1) = a + ppoly_coef(i1,2) = b + ppoly_coef(i1,3) = c + ppoly_coef(i1,4) = d + ppoly_coef(i1,5) = e end subroutine PQM_boundary_extrapolation_v1 diff --git a/src/ALE/polynomial_functions.F90 b/src/ALE/polynomial_functions.F90 index b0d5d135d5..0cc4eb0b71 100644 --- a/src/ALE/polynomial_functions.F90 +++ b/src/ALE/polynomial_functions.F90 @@ -21,62 +21,58 @@ module polynomial_functions contains ! ----------------------------------------------------------------------------- -! Pointwise evaluation of a polynomial -! ----------------------------------------------------------------------------- -real function evaluation_polynomial( coefficients, nb_coefficients, x ) +!> Pointwise evaluation of a polynomial at x +real function evaluation_polynomial( coeff, ncoef, x ) + real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial + integer, intent(in) :: ncoef !< The number of polynomial coefficients + real, intent(in) :: x !< The position at which to evaluate the polynomial ! ----------------------------------------------------------------------------- ! The polynomial is defined by the coefficients contained in the ! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ... -! where C refers to the array 'coefficients'. -! The number of coefficients is given by nb_coefficients and x +! where C refers to the array 'coeff'. +! The number of coefficients is given by ncoef and x ! is the coordinate where the polynomial is to be evaluated. ! ! The function returns the value of the polynomial at x. ! ----------------------------------------------------------------------------- ! Arguments - real, dimension(:), intent(in) :: coefficients - integer, intent(in) :: nb_coefficients - real, intent(in) :: x ! Local variables - integer :: k - real :: f ! value of polynomial at x + integer :: k + real :: f ! value of polynomial at x f = 0.0 - do k = 1,nb_coefficients - f = f + coefficients(k) * ( x**(k-1) ) + do k = 1,ncoef + f = f + coeff(k) * ( x**(k-1) ) end do evaluation_polynomial = f end function evaluation_polynomial -!> Calculates the first derivative of a polynomial with coefficients as above -!! evaluated at a point x -real function first_derivative_polynomial( coefficients, nb_coefficients, x ) +!> Calculates the first derivative of a polynomial evaluated at a point x +real function first_derivative_polynomial( coeff, ncoef, x ) + real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial + integer, intent(in) :: ncoef !< The number of polynomial coefficients + real, intent(in) :: x !< The position at which to evaluate the derivative ! ----------------------------------------------------------------------------- ! The polynomial is defined by the coefficients contained in the ! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ... -! where C refers to the array 'coefficients'. -! The number of coefficients is given by nb_coefficients and x +! where C refers to the array 'coeff'. +! The number of coefficients is given by ncoef and x ! is the coordinate where the polynomial's derivative is to be evaluated. ! -! The function returns the value of the polynomial at x. +! The function returns the first derivative of the polynomial at x. ! ----------------------------------------------------------------------------- - ! Arguments - real, dimension(:), intent(in) :: coefficients - integer, intent(in) :: nb_coefficients - real, intent(in) :: x - ! Local variables integer :: k real :: f ! value of polynomial at x f = 0.0 - do k = 2,nb_coefficients - f = f + REAL(k-1)*coefficients(k) * ( x**(k-2) ) + do k = 2,ncoef + f = f + REAL(k-1)*coeff(k) * ( x**(k-2) ) end do first_derivative_polynomial = f @@ -84,48 +80,45 @@ real function first_derivative_polynomial( coefficients, nb_coefficients, x ) end function first_derivative_polynomial ! ----------------------------------------------------------------------------- -! Exact integration of polynomial of degree n -! ----------------------------------------------------------------------------- -real function integration_polynomial( xi0, xi1, C, n ) +!> Exact integration of polynomial of degree npoly +real function integration_polynomial( xi0, xi1, Coeff, npoly ) + real, intent(in) :: xi0 !< The lower bound of the integral + real, intent(in) :: xi1 !< The lower bound of the integral + real, dimension(:), intent(in) :: Coeff !< The coefficients of the polynomial + integer, intent(in) :: npoly !< The degree of the polynomial ! ----------------------------------------------------------------------------- -! Exact integration of a polynomial of degree n over the interval [xi0,xi1]. -! The array of coefficients (C) must be of size n+1, where n is the degree of -! the polynomial to integrate. +! Exact integration of a polynomial of degree npoly over the interval [xi0,xi1]. +! The array of coefficients (Coeff) must be of size npoly+1. ! ----------------------------------------------------------------------------- - ! Arguments - real, intent(in) :: xi0, xi1 - real, dimension(:), intent(in) :: C - integer, intent(in) :: n - ! Local variables integer :: k real :: integral integral = 0.0 - do k = 1,(n+1) - integral = integral + C(k) * (xi1**k - xi0**k) / real(k) + do k = 1,npoly+1 + integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) end do ! !One non-answer-changing way of unrolling the above is: ! k=1 -! integral = integral + C(k) * (xi1**k - xi0**k) / real(k) -! if (n>=1) then +! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) +! if (npoly>=1) then ! k=2 -! integral = integral + C(k) * (xi1**k - xi0**k) / real(k) +! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) ! endif -! if (n>=2) then +! if (npoly>=2) then ! k=3 -! integral = integral + C(k) * (xi1**k - xi0**k) / real(k) +! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) ! endif -! if (n>=3) then +! if (npoly>=3) then ! k=4 -! integral = integral + C(k) * (xi1**k - xi0**k) / real(k) +! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) ! endif -! if (n>=4) then +! if (npoly>=4) then ! k=5 -! integral = integral + C(k) * (xi1**k - xi0**k) / real(k) +! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) ! endif ! integration_polynomial = integral diff --git a/src/ALE/regrid_edge_slopes.F90 b/src/ALE/regrid_edge_slopes.F90 index f8781aa937..e07f3c3bd5 100644 --- a/src/ALE/regrid_edge_slopes.F90 +++ b/src/ALE/regrid_edge_slopes.F90 @@ -32,6 +32,13 @@ module regrid_edge_slopes !------------------------------------------------------------------------------ !> Compute ih4 edge slopes (implicit third order accurate) subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) + 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) :: edge_slopes !< Returned edge slopes, with the + !! same units as u divided by the units of h. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ----------------------------------------------------------------------------- ! Compute edge slopes based on third-order implicit estimates. Note that ! the estimates are fourth-order accurate on uniform grids @@ -58,15 +65,6 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) ! boundary conditions close the system. ! ----------------------------------------------------------------------------- - ! Arguments - 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) :: edge_slopes !< Returned edge slopes, with the - !! same units as u divided by the units of h. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. - ! Local variables integer :: i, j ! loop indexes real :: h0, h1 ! cell widths @@ -188,6 +186,13 @@ end subroutine edge_slopes_implicit_h3 !------------------------------------------------------------------------------ !> Compute ih5 edge values (implicit fifth order accurate) subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) + 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) :: edge_slopes !< Returned edge slopes, with the + !! same units as u divided by the units of h. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ----------------------------------------------------------------------------- ! Fifth-order implicit estimates of edge values are based on a four-cell, ! three-edge stencil. A tridiagonal system is set up and is based on @@ -221,15 +226,6 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) ! on nonuniform meshes turned out to be intractable. ! ----------------------------------------------------------------------------- - ! Arguments - 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) :: edge_slopes !< Returned edge slopes, with the - !! same units as u divided by the units of h. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. - ! Local variables integer :: i, j, k ! loop indexes real :: h0, h1, h2, h3 ! cell widths diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index 1df7c6ec69..f9cdad794a 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -43,9 +43,15 @@ module regrid_edge_values contains !------------------------------------------------------------------------------ -! Bound edge values by neighboring cell averages -!------------------------------------------------------------------------------ -subroutine bound_edge_values( N, h, u, edge_values, h_neglect ) +!> Bound edge values by neighboring cell averages +subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) + 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) :: edge_val !< Potentially modified edge values, + !! with the same units as u. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ------------------------------------------------------------------------------ ! In this routine, we loop on all cells to bound their left and right ! edge values by the cell averages. That is, the left edge value must lie @@ -57,15 +63,6 @@ subroutine bound_edge_values( N, h, u, edge_values, h_neglect ) ! Therefore, boundary cells are treated as if they were local extrama. ! ------------------------------------------------------------------------------ - ! Arguments - 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) :: edge_values !< Potentially modified edge values, - !! 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 integer :: k0, k1, k2 @@ -111,8 +108,8 @@ subroutine bound_edge_values( N, h, u, edge_values, h_neglect ) u_c = u(k1) u_r = u(k2) - u0_l = edge_values(k,1) - u0_r = edge_values(k,2) + u0_l = edge_val(k,1) + u0_r = edge_val(k,2) 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 ) @@ -143,8 +140,8 @@ subroutine bound_edge_values( N, h, u, edge_values, h_neglect ) u0_r = max( min( u0_r, max(u_r, u_c) ), min(u_r, u_c) ) ! Store edge values - edge_values(k,1) = u0_l - edge_values(k,2) = u0_r + edge_val(k,1) = u0_l + edge_val(k,2) = u0_r end do ! loop on interior edges @@ -152,18 +149,16 @@ end subroutine bound_edge_values !------------------------------------------------------------------------------ -! Average discontinuous edge values (systematically) -!------------------------------------------------------------------------------ -subroutine average_discontinuous_edge_values( N, edge_values ) +!> Replace discontinuous collocated edge values with their average +subroutine average_discontinuous_edge_values( N, edge_val ) + integer, intent(in) :: N !< Number of cells + real, dimension(:,:), intent(inout) :: edge_val !< Edge values that may be modified; + !! the second index size is 2. ! ------------------------------------------------------------------------------ ! For each interior edge, check whether the edge values are discontinuous. -! If so, compute the average and replace the edge values by the average.! +! If so, compute the average and replace the edge values by the average. ! ------------------------------------------------------------------------------ - ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:,:), intent(inout) :: edge_values - ! Local variables integer :: k ! loop index real :: u0_minus ! left value at given edge @@ -174,15 +169,15 @@ subroutine average_discontinuous_edge_values( N, edge_values ) do k = 1,N-1 ! Edge value on the left of the edge - u0_minus = edge_values(k,2) + u0_minus = edge_val(k,2) ! Edge value on the right of the edge - u0_plus = edge_values(k+1,1) + u0_plus = edge_val(k+1,1) if ( u0_minus /= u0_plus ) then u0_avg = 0.5 * ( u0_minus + u0_plus ) - edge_values(k,2) = u0_avg - edge_values(k+1,1) = u0_avg + edge_val(k,2) = u0_avg + edge_val(k+1,1) = u0_avg end if end do ! end loop on interior edges @@ -191,19 +186,16 @@ end subroutine average_discontinuous_edge_values !------------------------------------------------------------------------------ -! Check discontinuous edge values and take average is not monotonic -!------------------------------------------------------------------------------ -subroutine check_discontinuous_edge_values( N, u, edge_values ) +!> Check discontinuous edge values and replace them with their average if not monotonic +subroutine check_discontinuous_edge_values( N, u, edge_val ) + integer, intent(in) :: N !< Number of cells + real, dimension(:), intent(in) :: u !< cell averages (size N) + real, dimension(:,:), intent(inout) :: edge_val !< Cell edge values with the same units as u. ! ------------------------------------------------------------------------------ ! For each interior edge, check whether the edge values are discontinuous. ! If so and if they are not monotonic, replace each edge value by their average. ! ------------------------------------------------------------------------------ - ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: u ! cell averages (size N) - real, dimension(:,:), intent(inout) :: edge_values - ! Local variables integer :: k ! loop index real :: u0_minus ! left value at given edge @@ -216,10 +208,10 @@ subroutine check_discontinuous_edge_values( N, u, edge_values ) do k = 1,N-1 ! Edge value on the left of the edge - u0_minus = edge_values(k,2) + u0_minus = edge_val(k,2) ! Edge value on the right of the edge - u0_plus = edge_values(k+1,1) + u0_plus = edge_val(k+1,1) ! Left cell average um_minus = u(k) @@ -230,8 +222,8 @@ subroutine check_discontinuous_edge_values( N, u, edge_values ) if ( (u0_plus - u0_minus)*(um_plus - um_minus) < 0.0 ) then u0_avg = 0.5 * ( u0_minus + u0_plus ) u0_avg = max( min( u0_avg, max(um_minus, um_plus) ), min(um_minus, um_plus) ) - edge_values(k,2) = u0_avg - edge_values(k+1,1) = u0_avg + edge_val(k,2) = u0_avg + edge_val(k+1,1) = u0_avg end if end do ! end loop on interior edges @@ -241,7 +233,14 @@ end subroutine check_discontinuous_edge_values !------------------------------------------------------------------------------ !> Compute h2 edge values (explicit second order accurate) -subroutine edge_values_explicit_h2( N, h, u, edge_values, h_neglect ) +subroutine edge_values_explicit_h2( N, h, u, edge_val, h_neglect ) + 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) :: edge_val !< Returned edge values, with the + !! same units as u; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ------------------------------------------------------------------------------ ! Compute edge values based on second-order explicit estimates. ! These estimates are based on a straight line spanning two cells and evaluated @@ -254,16 +253,7 @@ subroutine edge_values_explicit_h2( N, h, u, edge_values, h_neglect ) ! ! Boundary edge values are set to be equal to the boundary cell averages. ! ------------------------------------------------------------------------------ - - ! Arguments - 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) :: edge_values !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. - + ! Local variables integer :: k ! loop index real :: h0, h1 ! cell widths @@ -288,24 +278,31 @@ subroutine edge_values_explicit_h2( N, h, u, edge_values, h_neglect ) u1 = u(k) ! Compute left edge value - edge_values(k,1) = ( u0*h1 + u1*h0 ) / ( h0 + h1 ) + edge_val(k,1) = ( u0*h1 + u1*h0 ) / ( h0 + h1 ) ! Left edge value of the current cell is equal to right edge ! value of left cell - edge_values(k-1,2) = edge_values(k,1) + edge_val(k-1,2) = edge_val(k,1) end do ! end loop on interior cells ! Boundary edge values are simply equal to the boundary cell averages - edge_values(1,1) = u(1) - edge_values(N,2) = u(N) + edge_val(1,1) = u(1) + edge_val(N,2) = u(N) end subroutine edge_values_explicit_h2 !------------------------------------------------------------------------------ !> Compute h4 edge values (explicit fourth order accurate) -subroutine edge_values_explicit_h4( N, h, u, edge_values, h_neglect ) +subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect ) + 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) :: edge_val !< Returned edge values, with the + !! same units as u; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ----------------------------------------------------------------------------- ! Compute edge values based on fourth-order explicit estimates. ! These estimates are based on a cubic interpolant spanning four cells @@ -325,15 +322,6 @@ subroutine edge_values_explicit_h4( N, h, u, edge_values, h_neglect ) ! For this fourth-order scheme, at least four cells must exist. ! ----------------------------------------------------------------------------- - ! Arguments - 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) :: edge_values !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. - ! Local variables integer :: i, j real :: u0, u1, u2, u3 @@ -387,8 +375,8 @@ subroutine edge_values_explicit_h4( N, h, u, edge_values, h_neglect ) e = e / ( h0 + h1 + h2 + h3) - edge_values(i,1) = e - edge_values(i-1,2) = e + edge_val(i,1) = e + edge_val(i-1,2) = e #ifdef __DO_SAFETY_CHECKS__ if (e /= e) then @@ -422,14 +410,14 @@ subroutine edge_values_explicit_h4( N, h, u, edge_values, h_neglect ) call solve_linear_system( A, B, C, 4 ) ! First edge value - edge_values(1,1) = evaluation_polynomial( C, 4, x(1) ) + edge_val(1,1) = evaluation_polynomial( C, 4, x(1) ) ! Second edge value - edge_values(1,2) = evaluation_polynomial( C, 4, x(2) ) - edge_values(2,1) = edge_values(1,2) + edge_val(1,2) = evaluation_polynomial( C, 4, x(2) ) + edge_val(2,1) = edge_val(1,2) #ifdef __DO_SAFETY_CHECKS__ - if (edge_values(1,1) /= edge_values(1,1) .or. edge_values(1,2) /= edge_values(1,2)) then + if (edge_val(1,1) /= edge_val(1,1) .or. edge_val(1,2) /= edge_val(1,2)) then write(0,*) 'NaN in explicit_edge_h4 at k=',1 write(0,*) 'A=',A write(0,*) 'B=',B @@ -460,14 +448,14 @@ subroutine edge_values_explicit_h4( N, h, u, edge_values, h_neglect ) call solve_linear_system( A, B, C, 4 ) ! Last edge value - edge_values(N,2) = evaluation_polynomial( C, 4, x(5) ) + edge_val(N,2) = evaluation_polynomial( C, 4, x(5) ) ! Second to last edge value - edge_values(N,1) = evaluation_polynomial( C, 4, x(4) ) - edge_values(N-1,2) = edge_values(N,1) + edge_val(N,1) = evaluation_polynomial( C, 4, x(4) ) + edge_val(N-1,2) = edge_val(N,1) #ifdef __DO_SAFETY_CHECKS__ - if (edge_values(N,1) /= edge_values(N,1) .or. edge_values(N,2) /= edge_values(N,2)) then + if (edge_val(N,1) /= edge_val(N,1) .or. edge_val(N,2) /= edge_val(N,2)) then write(0,*) 'NaN in explicit_edge_h4 at k=',N write(0,*) 'A=' do i = 1,4 @@ -490,7 +478,14 @@ end subroutine edge_values_explicit_h4 !------------------------------------------------------------------------------ !> Compute ih4 edge values (implicit fourth order accurate) -subroutine edge_values_implicit_h4( N, h, u, edge_values, h_neglect ) +subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) + 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) :: edge_val !< Returned edge values, with the + !! same units as u; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ----------------------------------------------------------------------------- ! Compute edge values based on fourth-order implicit estimates. ! @@ -515,15 +510,6 @@ subroutine edge_values_implicit_h4( N, h, u, edge_values, h_neglect ) ! boundary conditions close the system. ! ----------------------------------------------------------------------------- - ! Arguments - 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) :: edge_values !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. - ! Local variables integer :: i, j ! loop indexes real :: h0, h1 ! cell widths @@ -627,18 +613,25 @@ subroutine edge_values_implicit_h4( N, h, u, edge_values, h_neglect ) call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) do i = 2,N - edge_values(i,1) = tri_x(i) - edge_values(i-1,2) = tri_x(i) + edge_val(i,1) = tri_x(i) + edge_val(i-1,2) = tri_x(i) end do - edge_values(1,1) = tri_x(1) - edge_values(N,2) = tri_x(N+1) + edge_val(1,1) = tri_x(1) + edge_val(N,2) = tri_x(N+1) end subroutine edge_values_implicit_h4 !------------------------------------------------------------------------------ !> Compute ih6 edge values (implicit sixth order accurate) -subroutine edge_values_implicit_h6( N, h, u, edge_values, h_neglect ) +subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) + 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) :: edge_val !< Returned edge values, with the + !! same units as u; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width + !! in the same units as h. ! ----------------------------------------------------------------------------- ! Sixth-order implicit estimates of edge values are based on a four-cell, ! three-edge stencil. A tridiagonal system is set up and is based on @@ -672,15 +665,6 @@ subroutine edge_values_implicit_h6( N, h, u, edge_values, h_neglect ) ! on nonuniform meshes turned out to be intractable. ! ----------------------------------------------------------------------------- - ! Arguments - 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) :: edge_values !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. - ! Local variables integer :: i, j, k ! loop indexes real :: h0, h1, h2, h3 ! cell widths @@ -1124,11 +1108,11 @@ subroutine edge_values_implicit_h6( N, h, u, edge_values, h_neglect ) call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) do i = 2,N - edge_values(i,1) = tri_x(i) - edge_values(i-1,2) = tri_x(i) + edge_val(i,1) = tri_x(i) + edge_val(i-1,2) = tri_x(i) end do - edge_values(1,1) = tri_x(1) - edge_values(N,2) = tri_x(N+1) + edge_val(1,1) = tri_x(1) + edge_val(N,2) = tri_x(N+1) end subroutine edge_values_implicit_h6 diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 6858e0cded..d9d2a19228 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -18,8 +18,7 @@ module regrid_interp implicit none ; private -type, public :: interp_CS_type - private +type, public :: interp_CS_type ; private !> The following parameter is only relevant when used with the target !! interface densities regridding scheme. It indicates which interpolation @@ -476,7 +475,9 @@ end function get_polynomial_coordinate !> Numeric value of interpolation_scheme corresponding to scheme name integer function interpolation_scheme(interp_scheme) - character(len=*), intent(in) :: interp_scheme !< Name of interpolation scheme + character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme + !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4", + !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5" select case ( uppercase(trim(interp_scheme)) ) case ("P1M_H2"); interpolation_scheme = INTERPOLATION_P1M_H2 @@ -494,18 +495,23 @@ integer function interpolation_scheme(interp_scheme) end select end function interpolation_scheme +!> Store the interpolation_scheme value in the interp_CS based on the input string. subroutine set_interp_scheme(CS, interp_scheme) - type(interp_CS_type), intent(inout) :: CS - character(len=*), intent(in) :: interp_scheme + type(interp_CS_type), intent(inout) :: CS !< A control structure for regrid_interp + character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme + !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4", + !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5" CS%interpolation_scheme = interpolation_scheme(interp_scheme) end subroutine set_interp_scheme -subroutine set_interp_extrap(CS, extrapolation) - type(interp_CS_type), intent(inout) :: CS - logical, intent(in) :: extrapolation +!> Store the boundary_extrapolation value in the interp_CS +subroutine set_interp_extrap(CS, extrap) + type(interp_CS_type), intent(inout) :: CS !< A control structure for regrid_interp + logical, intent(in) :: extrap !< Indicate whether high-order boundary + !! extrapolation should be used in boundary cells - CS%boundary_extrapolation = extrapolation + CS%boundary_extrapolation = extrap end subroutine set_interp_extrap end module regrid_interp diff --git a/src/ALE/regrid_solvers.F90 b/src/ALE/regrid_solvers.F90 index 1f15f97d1b..7e44039831 100644 --- a/src/ALE/regrid_solvers.F90 +++ b/src/ALE/regrid_solvers.F90 @@ -25,21 +25,18 @@ module regrid_solvers contains ! ----------------------------------------------------------------------------- -! Solve the linear system AX = B -! ----------------------------------------------------------------------------- +!> Solve the linear system AX = B by Gaussian elimination subroutine solve_linear_system( A, B, X, system_size ) + real, dimension(:,:), intent(inout) :: A !< The matrix being inverted + real, dimension(:), intent(inout) :: B !< system right-hand side + real, dimension(:), intent(inout) :: X !< solution vector + integer, intent(in) :: system_size !< The size of the system ! ----------------------------------------------------------------------------- ! This routine uses Gauss's algorithm to transform the system's original ! matrix into an upper triangular matrix. Back substitution yields the answer. ! The matrix A must be square and its size must be that of the vectors B and X. ! ----------------------------------------------------------------------------- - ! Arguments - real, dimension(:,:), intent(inout) :: A - real, dimension(:), intent(inout) :: B - real, dimension(:), intent(inout) :: X - integer :: system_size - ! Local variables integer :: i, j, k real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed @@ -127,18 +124,18 @@ end subroutine solve_linear_system ! ----------------------------------------------------------------------------- -! Solve the tridiagonal system AX = B -! ----------------------------------------------------------------------------- +!> Solve the tridiagonal system AX = B subroutine solve_tridiagonal_system( Al, Ad, Au, B, X, system_size ) + real, dimension(:), intent(inout) :: Ad !< Maxtix center diagonal + real, dimension(:), intent(inout) :: Al !< Matrix lower diagonal + real, dimension(:), intent(inout) :: Au !< Matrix upper diagonal + real, dimension(:), intent(inout) :: B !< system right-hand side + real, dimension(:), intent(inout) :: X !< solution vector + integer, intent(in) :: system_size !< The size of the system ! ----------------------------------------------------------------------------- ! This routine uses Thomas's algorithm to solve the tridiagonal system AX = B. ! (A is made up of lower, middle and upper diagonals) ! ----------------------------------------------------------------------------- - ! Arguments - real, dimension(:), intent(inout) :: Al, Ad, Au ! lo., mid. and up. diagonals - real, dimension(:), intent(inout) :: B ! system right-hand side - real, dimension(:), intent(inout) :: X ! solution vector - integer, intent(in) :: system_size ! Local variables integer :: k ! Loop index