Skip to content
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

Adding area correction based on 'epsilon bubble' for RRM (redo PR1022) #1581

Merged
merged 4 commits into from
Jul 2, 2017
+247 −39
Merged

Adding area correction based on 'epsilon bubble' for RRM (redo PR1022) #1581

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
53ecf17
Adding area correction based on 'epsilon bubble' for RRM (redo PR1022)
oksanaguba Jun 13, 2017
34eb159
Fixes issue that repro_sum_mod is not in CAM.
oksanaguba Jun 15, 2017
5000116
Fixing land IC file.
oksanaguba Jun 15, 2017
97f454d
cleaned up 'gp' code per questions left in comments
mt5555 Jul 1, 2017
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
8 changes: 4 additions & 4 deletions cime/config/acme/config_grids.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1451,10 +1451,10 @@
<nx>89147</nx>
<ny>1</ny>
<desc>1-deg with 1/4-deg over Cont. U.S. (version 1):</desc>
<file atm_mask="tx0.1v2">domain.lnd.conusx4v1_tx0.1v2.141022.nc</file>
<file ice_mask="tx0.1v2">domain.ocn.conusx4v1_tx0.1v2.141022.nc</file>
<file lnd_mask="tx0.1v2">domain.lnd.conusx4v1_tx0.1v2.141022.nc</file>
<file ocn_mask="tx0.1v2">domain.ocn.conusx4v1_tx0.1v2.141022.nc</file>
<file atm_mask="tx0.1v2">domain.lnd.conusx4v1_tx0.1v2.161129.nc</file>
<file ice_mask="tx0.1v2">domain.ocn.conusx4v1_tx0.1v2.161129.nc</file>
<file lnd_mask="tx0.1v2">domain.lnd.conusx4v1_tx0.1v2.161129.nc</file>
<file ocn_mask="tx0.1v2">domain.ocn.conusx4v1_tx0.1v2.161129.nc</file>
</domain>

<domain name="ne0np4_svalbard_x8v1_lowcon">
Expand Down
2 changes: 1 addition & 1 deletion components/cam/bld/namelist_files/namelist_defaults_cam.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1477,7 +1477,7 @@
<finidat hgrid="ne0np4_arm_x8v3_lowcon" ic_ymd="101" sim_year="2000">lnd/clm2/initdata/clmi.armx8v3.1850-01-01.nc</finidat>
<!-- for CLM4.0finidat hgrid="ne0np4_conus_x4v1_lowcon" ic_ymd="101" sim_year="2000">lnd/clm2/initdata/clmi.conusx4v1.2000-01-01_c141106.nc</finidat-->
<finidat hgrid="ne0np4_conus_x4v1_lowcon" ic_ymd="101" sim_year="1850">lnd/clm2/initdata_map/clmi.I1850CLM45.conusx4v1.74e105b.clm2.r.0021-01-01-00000.nc</finidat>
<finidat hgrid="ne0np4_conus_x4v1_lowcon" ic_ymd="101" sim_year="2000">lnd/clm2/initdata_map/clmi.ICRUCLM45.conusx4v1.74e105b.clm2.r.0021-01-01-00000.nc</finidat>
<finidat hgrid="ne0np4_conus_x4v1_lowcon" ic_ymd="101" sim_year="2000">lnd/clm2/initdata_map/clmi.ICRUCLM45.0021-01-01.conusx4v1.0ac6464_simyr2000_c161130.nc</finidat>
<finidat hgrid="ne0np4_svalbard_x8v1_lowcon" ic_ymd="101" sim_year="1850">lnd/clm2/initdata/clmi.svalbardx8v1.1850-01-01.nc</finidat>
<finidat hgrid="ne0np4_sooberingoa_x4x8v1_lowcon" ic_ymd="101" sim_year="1850">lnd/clm2/initdata/clmi.sooberingoax4x8v1.1850-01-01_c141110.nc</finidat>

Expand Down
193 changes: 186 additions & 7 deletions components/homme/src/share/cube_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
!--------------------------------------------------------------------------------
!
! 08/2016: O. Guba Modifying metric_atomic routine to add support for 'epsilon
! bubble' reference element map with GLL area = geometric area
!

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
Expand Down Expand Up @@ -78,6 +84,8 @@ module cube_mod
public :: CubeSetupEdgeIndex
public :: ref2sphere

public :: set_area_correction_map0, set_area_correction_map2

! public interface to REFERECE element map
#if HOMME_QUAD_PREC
interface ref2sphere
Expand Down Expand Up @@ -438,13 +446,37 @@ subroutine metric_atomic(elem,gll_points,alpha)
elem%dx_short = 1.0d0/(max_svd*0.5d0*dble(np-1)*rrearth*1000.0d0)
elem%dx_long = 1.0d0/(min_svd*0.5d0*dble(np-1)*rrearth*1000.0d0)

! optional noramlization:
elem%D = elem%D * sqrt(alpha)
elem%Dinv = elem%Dinv / sqrt(alpha)
elem%metdet = elem%metdet * alpha
elem%rmetdet = elem%rmetdet / alpha
elem%met = elem%met * alpha
elem%metinv = elem%metinv / alpha
! Area correction: Bring numerical area from integration weights to
! agreement with geometric area.
! Three different cases:
! (1) alpha == 1, this means that cube_init_atomic wasn't
! called with alpha parameter there will be no correction,
! (2) alpha <> 1 and cubed_sphere_map=0 and correction is so-called
! 'alpha-correction',
! (3) alpha <> 1 and cubed_sphere_map=2 and it is 'epsilon-bubble'
! correction.

if( cubed_sphere_map == 0 ) then
! alpha correction for cases (1) and (2).
elem%D = elem%D * sqrt(alpha)
elem%Dinv = elem%Dinv / sqrt(alpha)
elem%metdet = elem%metdet * alpha
elem%rmetdet = elem%rmetdet / alpha
elem%met = elem%met * alpha
elem%metinv = elem%metinv / alpha
elseif( cubed_sphere_map == 2 ) then
! eps bubble correction for case (3).
do j=2,np-1
do i=2,np-1
elem%D(i,j,:,:) = elem%D(i,j,:,:) * sqrt(alpha)
elem%Dinv(i,j,:,:) = elem%Dinv(i,j,:,:) / sqrt(alpha)
elem%metdet(i,j) = elem%metdet(i,j) * alpha
elem%rmetdet(i,j) = elem%rmetdet(i,j) / alpha
elem%met(i,j,:,:) = elem%met(i,j,:,:) * alpha
elem%metinv(i,j,:,:) = elem%metinv(i,j,:,:) / alpha
enddo
enddo
endif ! end of alpha/eps. bubble correction

end subroutine metric_atomic

Expand Down Expand Up @@ -2058,6 +2090,153 @@ function ref2sphere_elementlocal_q(q, corners) result(sphere)
end function


subroutine set_area_correction_map0(elem, nelemd, par, gp)
! Numerical area of the domain (sphere) is sum of integration weights. The sum
! is not exactly equal geometric area (4\pi R). It is required that
! numerical area = geometric area.
! This correction 'butters'
! the difference between numerical and geometrical areas evenly among DOFs. Then
! geometric areas of individual elements still do not equal their numerical areas,
! only whole domain's areas coinside.


#ifndef CAM
use repro_sum_mod, only: repro_sum
#else
use shr_reprosum_mod, only: repro_sum => shr_reprosum_calc
#endif

use quadrature_mod, only: quadrature_t
use element_mod, only: element_t
use parallel_mod, only: parallel_t
use kinds, only: iulog
use control_mod, only: topology

implicit none

type (element_t), pointer, intent(inout) :: elem(:)
type (parallel_t), intent(in) :: par
integer, intent(in) :: nelemd
type (quadrature_t), intent(in) :: gp

real(kind=real_kind), allocatable :: aratio(:,:)
real(kind=real_kind) :: area(1)
integer :: ie, i, j

allocate(aratio(nelemd,1))

if ( topology == "cube" ) then
area = 0.0d0
do ie=1,nelemd
aratio(ie,1) = sum(elem(ie)%mp(:,:)*elem(ie)%metdet(:,:))
enddo

call repro_sum(aratio, area, nelemd, nelemd, 1, commid=par%comm)
area(1) = 4*dd_pi/area(1) ! ratio correction

if (par%masterproc) &
write(iulog,'(a,f20.17)') " re-initializing cube elements: alpha area correction=",&
area(1)

do ie=1,nelemd
call cube_init_atomic(elem(ie),gp%points,area(1))
enddo
endif ! end of topology == 'cube'
deallocate(aratio)

end subroutine set_area_correction_map0


subroutine set_area_correction_map2(elem, nelemd, par, gp)
! Numerical area of the domain (sphere) is sum of integration weights. The sum
! is not exactly equal geometric area (4\pi R). It is required that
! numerical area = geometric area.
! The 'epsilon bubble' approach modifies inner weights in each element so that
! geometic and numerical areas of each element match.
#ifndef CAM
use repro_sum_mod, only: repro_sum
#else
use shr_reprosum_mod, only: repro_sum => shr_reprosum_calc
#endif

use quadrature_mod, only: quadrature_t
use element_mod, only: element_t
use parallel_mod, only: parallel_t
use kinds, only: iulog
use control_mod, only: topology

implicit none

type (element_t), pointer, intent(inout) :: elem(:)
type (parallel_t), intent(in) :: par
integer, intent(in) :: nelemd
type (quadrature_t), intent(in) :: gp

real(kind=real_kind), allocatable :: aratio(:,:)
real(kind=real_kind) :: area(1), area_sphere, area_num, area_dummy,sum_w, delta
integer :: ie, i, j
real(kind=real_kind) :: tol_local = 1e-15

allocate(aratio(nelemd,1))

if ( topology == "cube" ) then
do ie=1,nelemd
! Obtain area of element = sum of areas of 2 triangles.
call sphere_tri_area(elem(ie)%corners3D(1), elem(ie)%corners3D(2),&
elem(ie)%corners3D(3), area_sphere)
call sphere_tri_area(elem(ie)%corners3D(1), elem(ie)%corners3D(3),&
elem(ie)%corners3D(4), area_dummy)
! Store element's area in area_sphere.
area_sphere = area_sphere + area_dummy

! Compute 'numerical area' of the element as sum of integration
! weights.
area_num = 0.0d0
do j = 1,np
do i = 1,np
area_num = area_num + gp%weights(i)*gp%weights(j)*elem(ie)%metdet(i,j)
enddo
enddo

! Compute sum of inner integration weights for correction.
sum_w = 0.0d0 ! or sum_w = sum(elem(ie)%mp(2:np-1,2:np-1)*elem(ie)%metdet(2:np-1,2:np-1))
do j = 2, np-1
do i = 2, np-1
sum_w = sum_w + gp%weights(i)*gp%weights(j)*elem(ie)%metdet(i,j)
enddo
enddo
! Which tol is to use here?
if ( sum_w > tol_local ) then
delta = (area_sphere - area_num)/sum_w
call cube_init_atomic(elem(ie),gp%points,1.0d0 + delta)
else
! Abort since the denominator in correction is too small.
call abortmp('Cube_mod, set_area_correction_map2() : Area correction based on eps. bubble &
cannot be done, sum_w is too small.')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Getting this with the IBM compiler:

components/homme/src/share/cube_mod.F90", line 2215.1: 1515-021 (E) Syntax error: Token " & " is expected.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is only about this string then? Even for things this small, there should be PR? I will submit a PR for this (and another thing) tonight.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't have means to test this with IBm, is this an issue because the line is too long or because IBM does not like
print('........... &
..........') ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The ampersand at the beginning of the second/continuation line of a string appears to be mandatory for the IBM compiler:

                 call abortmp('Cube_mod, set_area_correction_map2() : Area correction based on eps. bubble &
-                              cannot be done, sum_w is too small.')
+                             &cannot be done, sum_w is too small.')
              endif

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or you break the string into two and not have the ampersand inside of a string.

endif
enddo ! loop over elements
! code for verification.
area = 0.0d0
do ie = 1,nelemd
aratio(ie,1) = 0.0
do j = 1,np
do i = 1,np
aratio(ie,1) = aratio(ie,1) + gp%weights(i)*gp%weights(j)*elem(ie)%metdet(i,j)
enddo
enddo
enddo

call repro_sum(aratio, area, nelemd, nelemd, 1, commid=par%comm)
if (par%masterproc) &
write(iulog,'(a,f20.17)') "Epsilon bubble correction: Corrected area - 4\pi ",area(1) - 4.0d0*dd_pi

endif ! end of topology == 'cube'
deallocate(aratio)

end subroutine set_area_correction_map2




end module cube_mod

Expand Down
45 changes: 23 additions & 22 deletions components/homme/src/share/prim_driver_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
! ------------------------------------------------------------------------------------------------
! prim_driver_mod:
!
! 08/2016: O. Guba Inserting code for "espilon bubble" reference element map
!
!
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
Expand Down Expand Up @@ -32,7 +38,6 @@ module prim_driver_mod
public :: smooth_topo_datasets

type (cg_t), allocatable :: cg(:) ! conjugate gradient struct (nthreads)
type (quadrature_t) :: gp ! element GLL points

#ifndef CAM
type (ColumnModel_t), allocatable :: cm(:) ! (nthreads)
Expand All @@ -48,7 +53,7 @@ subroutine prim_init1(elem, par, dom_mt, Tl)
use thread_mod, only : nthreads, nThreadsHoriz
! --------------------------------
use control_mod, only : runtype, restartfreq, integration, topology, &
partmethod, use_semi_lagrange_transport
partmethod, use_semi_lagrange_transport, cubed_sphere_map
! --------------------------------
use prim_state_mod, only : prim_printstate_init
! --------------------------------
Expand All @@ -60,11 +65,12 @@ subroutine prim_init1(elem, par, dom_mt, Tl)
! --------------------------------
use mass_matrix_mod, only : mass_matrix
! --------------------------------
use cube_mod, only : cubeedgecount , cubeelemcount, cubetopology
use cube_mod, only : cubeedgecount , cubeelemcount, cubetopology, cube_init_atomic, &
set_corner_coordinates, assign_node_numbers_to_elem, &
set_area_correction_map0, set_area_correction_map2
! --------------------------------
use mesh_mod, only : MeshSetCoordinates, MeshUseMeshFile, MeshCubeTopology, &
MeshCubeElemCount, MeshCubeEdgeCount
use cube_mod, only : cube_init_atomic, set_corner_coordinates, assign_node_numbers_to_elem
! --------------------------------
use metagraph_mod, only : metavertex_t, metaedge_t, localelemcount, initmetagraph, printmetavertex
! --------------------------------
Expand Down Expand Up @@ -132,16 +138,14 @@ subroutine prim_init1(elem, par, dom_mt, Tl)
integer :: err, ierr, l, j
logical, parameter :: Debug = .FALSE.

real(kind=real_kind), allocatable :: aratio(:,:)
real(kind=real_kind) :: area(1),xtmp

integer :: i
integer,allocatable :: TailPartition(:)
integer,allocatable :: HeadPartition(:)

integer total_nelem
real(kind=real_kind) :: approx_elements_per_task
integer :: n_domains
type (quadrature_t) :: gp ! element GLL points

#ifndef CAM
logical :: repro_sum_use_ddpdd, repro_sum_recompute
Expand Down Expand Up @@ -358,23 +362,20 @@ subroutine prim_init1(elem, par, dom_mt, Tl)
! =================================================================
if(par%masterproc) write(iulog,*) 'running mass_matrix'
call mass_matrix(par,elem)
allocate(aratio(nelemd,1))

if (topology=="cube") then
area = 0
do ie=1,nelemd
aratio(ie,1) = sum(elem(ie)%mp(:,:)*elem(ie)%metdet(:,:))
enddo
call repro_sum(aratio, area, nelemd, nelemd, 1, commid=par%comm)
area(1) = 4*dd_pi/area(1) ! ratio correction
deallocate(aratio)
if (par%masterproc) &
write(iulog,'(a,f20.17)') " re-initializing cube elements: area correction=",area(1)
! global area correction (default for cubed-sphere meshes)
if( cubed_sphere_map == 0 ) then
call set_area_correction_map0(elem, nelemd, par, gp)
endif

! Epsilon bubble correction (default for RRM meshes).
if(( cubed_sphere_map == 2 ).AND.( np > 2 )) then
call set_area_correction_map2(elem, nelemd, par, gp)
endif

deallocate(gp%points)
deallocate(gp%weights)

do ie=1,nelemd
call cube_init_atomic(elem(ie),gp%points,area(1))
enddo
end if

if(par%masterproc) write(iulog,*) 're-running mass_matrix'
call mass_matrix(par,elem)
Expand Down
Loading