From c41ea1963f69aab1682a51b72215bec91b01d6ab Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Sat, 16 Jan 2016 14:05:06 -0500 Subject: [PATCH 1/5] Added re-factor of remapByProjection() - New function remap_via_sub_cells() returns to the idea of remapping by projection implemented via sub-cell integrals with a correction step to the sub-cell integrals that ensures accurate conservation. - Also added a function for evaluating the average of a single cell's reconstruction over an interval within the cell. --- src/ALE/MOM_remapping.F90 | 222 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 222 insertions(+) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index a4ffeb3083..273887f9ff 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -426,6 +426,228 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) end subroutine remapping_core +!> Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating +!! the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the +!! appropriate integrals into the h1*u1 values. +subroutine remap_by_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h1, method, u1 ) + integer, intent(in) :: n0 !< Number of cells in source grid + real, intent(in) :: h0(:) !< Source grid widths (size n0) + real, intent(in) :: u0(:) !< Source cell averages (size n0) + real, intent(in) :: ppoly0_E(:,:) !< Edge value of polynomial + real, intent(in) :: ppoly0_coefficients(:,:) !< Coefficients of polynomial + integer, intent(in) :: n1 !< Number of cells in target grid + real, intent(in) :: h1(:) !< Target grid widths (size n1) + integer, intent(in) :: method !< Remapping scheme to use + real, intent(out) :: u1(:) !< Target cell averages (size n1) + ! Local variables + integer :: i_sub ! Index of sub-cell + integer :: i0 ! Index into h0(1:n0), source column + integer :: i1 ! Index into h1(1:n1), target column + integer :: i_start0 ! Used to record which sub-cells map to source cells + integer :: i_start1 ! Used to record which sub-cells map to target cells + integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell + real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell + real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell + real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell + real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell + integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell + integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell + integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell + integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell + integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell + real :: xa, xb ! Non-dimensional position within a source cell (0..1) + real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells + real :: dh ! The width of the sub-cell + real :: duh ! The total amount of accumulated stuff (u*h) + + ! Initialize algorithm + h0_supply = h0(1) + h1_supply = h1(1) + i0 = 1 ; i1 = 1 + xa = 0. + + ! First sub-cell is always vanished + h_sub(1) = 0. + isrc_start(1) = 1 + isrc_end(1) = 1 + isrc_max(1) = 1 + uh_sub(1) = 0. + u_sub(1) = u0(1) + + i_start0 = 1 + i_start1 = 1 + i_max = 1 + dh_max = 0. + + ! Loop over each sub-cell + do i_sub = 2, n0+n1+1 + + ! This is the width of the sub-cell, determined by which ever column has the least + ! supply available to consume. + dh = min(h0_supply, h1_supply) + + ! Record the source index (i0) that this sub-cell integral belongs to. This + ! is needed to index the reconstruction coefficients for the source cell + ! used in the integrals of the sub-cell width. + h_sub(i_sub) = dh + + ! For recording the largest sub-cell within a source cell. + if (dh >= dh_max) then + i_max = i_sub + dh_max = dh + endif + + ! Evaluate average and integral for sub-cell i_sub. + ! Integral is over distance dh but expressed in terms of non-dimensional + ! positions with source cell from xa to xb (0 <= xa <= xb <= 1). + xb = xa + dh / h0(i0) ! This expression should ideally yield xa <= xb <= 1.0 + !xb = min(xb, 1.0) ! BUT the last inequality could be violated due to roundoff - AJA + u_sub(i_sub) = average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0, xa, xb, dh) + uh_sub(i_sub) = dh * u_sub(i_sub) + xa = xb ! Next integral will start at end of last + + ! Which ever column (source or target) has the least width left to consume determined + ! the width, dh, of sub-cell i_sub in the expression for dh above. + if (h0_supply <= h1_supply) then + ! h0_supply is smaller than h1_supply) so we consume h0_supply and increment the + ! source cell index. + h1_supply = h1_supply - dh ! Although this is a difference the result will + ! be non-negative because of the conditional. + ! Record the sub-cell start/end index that span the source cell i0. + isrc_start(i0) = i_start0 + isrc_end(i0) = i_sub + i_start0 = i_sub + 1 + ! Record the sub-cell that is the largest fraction of the source cell. + isrc_max(i0) = i_max + i_max = i_sub + 1 + dh_max = 0. + ! Move the source index. + if (i0 < n0) then + i0 = i0 + 1 + h0_supply = h0(i0) + xa = 0. ! Next integral will start at top of the next source cell. + else + h0_supply = 0. + endif + else + ! h1_supply is smaller than h0_supply) so we consume h1_supply and increment the + ! target cell index. + h0_supply = h0_supply - dh ! Although this is a difference the result will + ! be non-negative because of the conditional. + ! Record the sub-cell start/end index that span the target cell i1. + itgt_start(i1) = i_start1 + itgt_end(i1) = i_sub + i_start1 = i_sub + 1 + ! Move the target index. + if (i1 < n1) then + i1 = i1 + 1 + h1_supply = h1(i1) + else + h1_supply = 0. + endif + endif + + enddo + + ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within + ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals + ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (@Hallberg-NOAA). + do i0 = 1, n0 + i_max = isrc_max(i0) + dh_max = h_sub(i_max) + if (dh_max > 0.) then + ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell. + duh = 0. + do i_sub = isrc_start(i0), isrc_end(i0) + if (i_sub /= i_max) duh = duh + uh_sub(i_sub) + enddo + uh_sub(i_max) = u0(i0)*h0(i0) - duh + u_sub(i_max) = uh_sub(i_max) / dh_max + endif + enddo + + ! Loop over each target cell summing the integrals from sub-cells within the target cell. + do i1 = 1, n1 + if (h1(i1) > 0.) then + duh = 0. + do i_sub = itgt_start(i1), itgt_end(i1) + duh = duh + uh_sub(i_sub) + enddo + u1(i1) = duh / h1(i1) + else + u1(i1) = u_sub(itgt_start(i1)) + endif + enddo + +end subroutine remap_by_via_sub_cells + +!> Returns the average value of a reconstruction within a single source cell, i0, +!! between the non-dimensional positions xa and xb (xa<=xb) with dimensional +!! separation dh. +real function average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0, xa, xb, dh) + integer, intent(in) :: n0 !< Number of cells in source grid + real, intent(in) :: ppoly0_E(:,:) !< Edge value of polynomial + real, intent(in) :: ppoly0_coefficients(:,:) !< Coefficients of polynomial + integer, intent(in) :: method !< Remapping scheme to use + integer, intent(in) :: i0 !< Source cell index + real, intent(in) :: xa !< Non-dimensional start position within source cell + real, intent(in) :: xb !< Non-dimensional end position within source cell + real, intent(in) :: dh !< Width of interval corresponding to xa..xb + ! Local variables + real :: u_ave, xa_2, xb_2, xa2pxb2, xapxb + real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials + + if (dh > 0.) then + select case ( method ) + case ( INTEGRATION_PCM ) + u_ave = ppoly0_coefficients(i0,1) + case ( INTEGRATION_PLM ) + u_ave = ( & + ppoly0_coefficients(i0,1) & + + ppoly0_coefficients(i0,2) * 0.5 * ( xb + xa ) ) + case ( INTEGRATION_PPM ) + u_ave = ( & + ppoly0_coefficients(i0,1) & + + ( ppoly0_coefficients(i0,2) * 0.5 * ( xb + xa ) & + + ppoly0_coefficients(i0,3) * r_3 * ( ( xb*xb + xa*xa ) + xa*xb ) ) ) + case ( INTEGRATION_PQM ) + xa_2 = xa*xa + xb_2 = xb*xb + xa2pxb2 = xa_2 + xb_2 + xapxb = xa + xb + u_ave = ( & + ppoly0_coefficients(i0,1) & + + ( ppoly0_coefficients(i0,2) * 0.5 * ( xapxb ) & + + ( ppoly0_coefficients(i0,3) * r_3 * ( xa2pxb2 + xa*xb ) & + + ( ppoly0_coefficients(i0,4) * 0.25* ( xa2pxb2 * xapxb ) & + + ppoly0_coefficients(i0,5) * 0.2 * ( ( xb*xb_2 + xa*xa_2 ) * xapxb + xa_2*xb_2 ) ) ) ) ) + case default + call MOM_error( FATAL,'The selected integration method is invalid' ) + end select + else ! dh == 0. + select case ( method ) + case ( INTEGRATION_PCM ) + u_ave = ppoly0_coefficients(i0,1) + case ( INTEGRATION_PLM ) + u_ave = ppoly0_coefficients(i0,1) & + + xa * ppoly0_coefficients(i0,2) + case ( INTEGRATION_PPM ) + u_ave = ppoly0_coefficients(i0,1) & + + xa * ( ppoly0_coefficients(i0,2) & + + xa * ppoly0_coefficients(i0,3) ) + case ( INTEGRATION_PQM ) + u_ave = ppoly0_coefficients(i0,1) & + + xa * ( ppoly0_coefficients(i0,2) & + + xa * ( ppoly0_coefficients(i0,3) & + + xa * ( ppoly0_coefficients(i0,4) & + + xa * ppoly0_coefficients(i0,5) ) ) ) + case default + call MOM_error( FATAL,'The selected integration method is invalid' ) + end select + endif + average_value_ppoly = u_ave + +end function average_value_ppoly !> Remaps column of values u0 on grid h0 to grid h1 by integrating !! over the projection of each h1 cell onto the h0 grid. From 9e60ae0f7e472c2064bd977dc8800e89d3b115d6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 18 Jan 2016 10:05:42 -0500 Subject: [PATCH 2/5] Fixed typo in s/r name remap_via_sub_cells - remap_by_via_sub_cells() is now remap_via_sub_cells() --- src/ALE/MOM_remapping.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 273887f9ff..7a2d4d14e6 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -429,7 +429,7 @@ end subroutine remapping_core !> Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating !! the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the !! appropriate integrals into the h1*u1 values. -subroutine remap_by_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h1, method, u1 ) +subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h1, method, u1 ) integer, intent(in) :: n0 !< Number of cells in source grid real, intent(in) :: h0(:) !< Source grid widths (size n0) real, intent(in) :: u0(:) !< Source cell averages (size n0) @@ -579,7 +579,7 @@ subroutine remap_by_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1 endif enddo -end subroutine remap_by_via_sub_cells +end subroutine remap_via_sub_cells !> Returns the average value of a reconstruction within a single source cell, i0, !! between the non-dimensional positions xa and xb (xa<=xb) with dimensional From 352dc31e5acba61480bf211acebbf75d3871876b Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 18 Jan 2016 10:38:00 -0500 Subject: [PATCH 3/5] Fixed target association for last sub-cell - On exiting the cycling over sub-cells, the remap_via_sub_cells() routine did not have an target association for the very last sub-cell. --- src/ALE/MOM_remapping.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 7a2d4d14e6..7762ce04df 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -549,6 +549,10 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h enddo + ! Last sub-cell needs to be associated with a target + itgt_start(n1) = i_start1 + itgt_end(n1) = n0+n1+1 + ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (@Hallberg-NOAA). From 4129446a85be00393a8e155d904082b326561015 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 18 Jan 2016 16:01:26 -0500 Subject: [PATCH 4/5] Debugged remapping to thicker target column - Separated integral calculations into their own loop - will later help create a vectorized version of remapping. - Expressed non-dimensional positions in terms of effective thickness so that positions are always bounded by source cell. - Removed unused argument, dh, to average_value_ppoly(). --- src/ALE/MOM_remapping.F90 | 83 ++++++++++++++++++++++++++------------- 1 file changed, 56 insertions(+), 27 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 7762ce04df..dd4217070f 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -450,45 +450,52 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell + integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell + real, dimension(n0) :: h0_eff ! Effective thickness of source cells integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell real :: xa, xb ! Non-dimensional position within a source cell (0..1) real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells real :: dh ! The width of the sub-cell real :: duh ! The total amount of accumulated stuff (u*h) + real :: dh0_eff ! Running sum of source cell thickness + real, parameter :: h_very_large = 1.E30 ! A large thickness, larger than will ever be encountered ! Initialize algorithm h0_supply = h0(1) h1_supply = h1(1) i0 = 1 ; i1 = 1 - xa = 0. + i_start0 = 1 ; i_start1 = 1 + i_max = 1 + dh_max = 0. + dh0_eff = 0. ! First sub-cell is always vanished h_sub(1) = 0. isrc_start(1) = 1 isrc_end(1) = 1 isrc_max(1) = 1 - uh_sub(1) = 0. - u_sub(1) = u0(1) - - i_start0 = 1 - i_start1 = 1 - i_max = 1 - dh_max = 0. + isub_src(1) = 1 - ! Loop over each sub-cell + ! Loop over each sub-cell to calculate intersections with source and target grids do i_sub = 2, n0+n1+1 ! This is the width of the sub-cell, determined by which ever column has the least ! supply available to consume. dh = min(h0_supply, h1_supply) + ! This is the running sum of the source cell thickness. After summing over each + ! sub-cell, the sum of sub-cell thickness might differ from the original source + ! cell thickness due to round off. + dh0_eff = dh0_eff + min(dh, h0_supply) + ! Record the source index (i0) that this sub-cell integral belongs to. This ! is needed to index the reconstruction coefficients for the source cell ! used in the integrals of the sub-cell width. + isub_src(i_sub) = i0 h_sub(i_sub) = dh ! For recording the largest sub-cell within a source cell. @@ -497,15 +504,6 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h dh_max = dh endif - ! Evaluate average and integral for sub-cell i_sub. - ! Integral is over distance dh but expressed in terms of non-dimensional - ! positions with source cell from xa to xb (0 <= xa <= xb <= 1). - xb = xa + dh / h0(i0) ! This expression should ideally yield xa <= xb <= 1.0 - !xb = min(xb, 1.0) ! BUT the last inequality could be violated due to roundoff - AJA - u_sub(i_sub) = average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0, xa, xb, dh) - uh_sub(i_sub) = dh * u_sub(i_sub) - xa = xb ! Next integral will start at end of last - ! Which ever column (source or target) has the least width left to consume determined ! the width, dh, of sub-cell i_sub in the expression for dh above. if (h0_supply <= h1_supply) then @@ -521,13 +519,15 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h isrc_max(i0) = i_max i_max = i_sub + 1 dh_max = 0. + ! Record the source cell thickness found by summing the sub-cell thicknesses. + h0_eff(i0) = dh0_eff ! Move the source index. if (i0 < n0) then i0 = i0 + 1 h0_supply = h0(i0) - xa = 0. ! Next integral will start at top of the next source cell. + dh0_eff = 0. else - h0_supply = 0. + h0_supply = h_very_large endif else ! h1_supply is smaller than h0_supply) so we consume h1_supply and increment the @@ -543,15 +543,45 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h i1 = i1 + 1 h1_supply = h1(i1) else - h1_supply = 0. + h1_supply = h_very_large endif endif enddo - ! Last sub-cell needs to be associated with a target - itgt_start(n1) = i_start1 - itgt_end(n1) = n0+n1+1 + ! Loop over each sub-cell to calculate average/integral values within each sub-cell. + xa = 0. + dh0_eff = 0. + uh_sub(1) = 0. + u_sub(1) = ppoly0_E(1,1) + do i_sub = 2, n0+n1 + + ! Sub-cell thickness from loop above + dh = h_sub(i_sub) + + ! Source cell + i0 = isub_src(i_sub) + + ! Evaluate average and integral for sub-cell i_sub. + ! Integral is over distance dh but expressed in terms of non-dimensional + ! positions with source cell from xa to xb (0 <= xa <= xb <= 1). + dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell + xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0 + xb = min(1., xb) ! This is only needed when the total target column is wider than the source column + u_sub(i_sub) = average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0, xa, xb) + uh_sub(i_sub) = dh * u_sub(i_sub) + + if (isub_src(i_sub+1) /= i0) then + ! If the next sub-cell is in a different source cell, reset the position counters + dh0_eff = 0. + xa = 0. + else + xa = xb ! Next integral will start at end of last + endif + + enddo + u_sub(n0+n1+1) = ppoly0_E(n0,2) ! This value is only needed when total target column + uh_sub(n0+n1+1) = ppoly0_E(n0,2) * h_sub(n0+n1+1) ! is wider than the source column ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals @@ -588,7 +618,7 @@ end subroutine remap_via_sub_cells !> Returns the average value of a reconstruction within a single source cell, i0, !! between the non-dimensional positions xa and xb (xa<=xb) with dimensional !! separation dh. -real function average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0, xa, xb, dh) +real function average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0, xa, xb) integer, intent(in) :: n0 !< Number of cells in source grid real, intent(in) :: ppoly0_E(:,:) !< Edge value of polynomial real, intent(in) :: ppoly0_coefficients(:,:) !< Coefficients of polynomial @@ -596,12 +626,11 @@ real function average_value_ppoly( n0, ppoly0_E, ppoly0_coefficients, method, i0 integer, intent(in) :: i0 !< Source cell index real, intent(in) :: xa !< Non-dimensional start position within source cell real, intent(in) :: xb !< Non-dimensional end position within source cell - real, intent(in) :: dh !< Width of interval corresponding to xa..xb ! Local variables real :: u_ave, xa_2, xb_2, xa2pxb2, xapxb real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials - if (dh > 0.) then + if (xb > xa) then select case ( method ) case ( INTEGRATION_PCM ) u_ave = ppoly0_coefficients(i0,1) From 6ae81b25e39fd005a4d1a927524033137eda102b Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 19 Jan 2016 09:32:43 -0500 Subject: [PATCH 5/5] Updated unit tests for MOM_remapping - remappingUnitTests now tests remap_via_sub_cells() --- src/ALE/MOM_remapping.F90 | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index dd4217070f..ff296871e6 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -1311,7 +1311,7 @@ subroutine end_remapping(CS) end subroutine end_remapping -!> Runs unit tests on remapping functions. +!> Runs unit tests on remapping functions. !! Should only be called from a single/root thread !! Returns True if a test fails, otherwise False logical function remappingUnitTests() @@ -1319,7 +1319,7 @@ logical function remappingUnitTests() real :: h0(n0), x0(n0+1), u0(n0) real :: h1(n1), x1(n1+1), u1(n1), hn1(n1), dx1(n1+1) real :: h2(n2), x2(n2+1), u2(n2), hn2(n2), dx2(n2+1) - data u0 /3., 1., -1., -3./ ! Linear profile, 4 at surface to -4 at bottom + data u0 /9., 3., -3., -9./ ! Linear profile, 4 at surface to -4 at bottom data h0 /4*0.75/ ! 4 uniform layers with total depth of 3 data h1 /3*1./ ! 3 uniform layers with total depth of 3 data h2 /6*0.5/ ! 6 uniform layers with total depth of 3 @@ -1356,7 +1356,7 @@ logical function remappingUnitTests() call dzFromH1H2( n0, h0, n1, h1, dx1 ) call remapping_core( CS, n0, h0, u0, n1, dx1, u1 ) do i=1,n1 - err=u1(i)-(8./3.)*(0.5*real(1+n1)-real(i)) + err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>real(n1-1)*epsilon(err)) thisTest = .true. enddo write(*,*) 'h1 (by projection)' @@ -1380,7 +1380,7 @@ logical function remappingUnitTests() call remapByProjection( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n1, h1, INTEGRATION_PPM, u1 ) do i=1,n1 - err=u1(i)-8./3.*(0.5*real(1+n1)-real(i)) + err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo if (thisTest) write(*,*) 'remappingUnitTests: Failed remapByProjection()' @@ -1395,7 +1395,7 @@ logical function remappingUnitTests() call dumpGrid(n1,h1,x1,u1) hn1=hn1-h1 do i=1,n1 - err=u1(i)-8./3.*(0.5*real(1+n1)-real(i)) + err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo if (thisTest) write(*,*) 'remappingUnitTests: Failed remapByDeltaZ() 1' @@ -1414,12 +1414,33 @@ logical function remappingUnitTests() call dumpGrid(n2,hn2,x2,u2) do i=1,n2 - err=u2(i)-8./6.*(0.5*real(1+n2)-real(i)) + err=u2(i)-8./2.*(0.5*real(1+n2)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo if (thisTest) write(*,*) 'remappingUnitTests: Failed remapByDeltaZ() 2' remappingUnitTests = remappingUnitTests .or. thisTest + write(*,*) 'Via sub-cells' + thisTest = .false. + call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & + n2, h2, INTEGRATION_PPM, u2 ) + call dumpGrid(n2,h2,x2,u2) + + do i=1,n2 + err=u2(i)-8./2.*(0.5*real(1+n2)-real(i)) + if (abs(err)>2.*epsilon(err)) thisTest = .true. + enddo + if (thisTest) write(*,*) 'remappingUnitTests: Failed remap_via_sub_cells() 2' + remappingUnitTests = remappingUnitTests .or. thisTest + + call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & + 6, (/.125,.125,.125,.125,.125,.125/), INTEGRATION_PPM, u2 ) + call dumpGrid(6,h2,x2,u2) + + call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & + 3, (/2.25,1.5,1./), INTEGRATION_PPM, u2 ) + call dumpGrid(3,h2,x2,u2) + deallocate(ppoly0_E, ppoly0_S, ppoly0_coefficients) write(*,*) '=========================================================='