Skip to content

Alternative remapping s/r #251

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Hallberg-NOAA merged 5 commits into from
Jan 19, 2016
Merged

Alternative remapping s/r #251

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c41ea19
Added re-factor of remapByProjection()
adcroft Jan 16, 2016
9e60ae0
Fixed typo in s/r name remap_via_sub_cells
adcroft Jan 18, 2016
352dc31
Fixed target association for last sub-cell
adcroft Jan 18, 2016
4129446
Debugged remapping to thicker target column
adcroft Jan 18, 2016
6ae81b2
Updated unit tests for MOM_remapping
adcroft Jan 19, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
288 changes: 282 additions & 6 deletions src/ALE/MOM_remapping.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -426,6 +426,261 @@ 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_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+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
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
isub_src(1) = 1

! 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.
if (dh >= dh_max) then
i_max = i_sub
dh_max = dh
endif

! 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.
! 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)
dh0_eff = 0.
else
h0_supply = h_very_large
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 = h_very_large
endif
endif

enddo

! 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
! 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_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)
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
! 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 (xb > xa) 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.
Expand Down Expand Up @@ -1056,15 +1311,15 @@ 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()
integer, parameter :: n0 = 4, n1 = 3, n2 = 6
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
Expand Down Expand Up @@ -1101,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)'
Expand All @@ -1125,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()'
Expand All @@ -1140,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'
Expand All @@ -1159,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(*,*) '=========================================================='
Expand Down