Skip to content

Commit 8364443

Browse files
OlgaSergienkoOlgaSergienko
committed
Adding metric terms to account for non-rectangular shape of grid cells
This PR corrects an implementation of a finite-element solver for ice-sheet/shelf velocity. It computes a Jacobian weight J_q in shape functions for quadrature points and applies it to integrals over cells. This correction is significant when the shape of the cell is not rectangular. Changes made: - nitialize_ice_shelf_dyn: allocates a new variable `Jac` - bilinear_shape_fn_grid: computes a Jacobian determinant (`Jac`) at each quadrature point - CG_action: uses `jac_wt= Jac * IareaT` for computing the ice deformation and basal traction terms - matrix_diagonal: uses `jac_wt= Jac * IareaT` for computing diagonal terms - ice_shelf_dyn_end: deallocates `Jac`
1 parent ebed2b4 commit 8364443

1 file changed

Lines changed: 28 additions & 11 deletions

File tree

‎src/ice_shelf/MOM_ice_shelf_dynamics.F90‎

Lines changed: 28 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,10 @@ module MOM_ice_shelf_dynamics
126126
!! 4 quadrature points surrounding the cell vertices [L-1 ~> m-1].
127127
real, pointer, dimension(:,:,:) :: PhiC => NULL() !< The gradients of bilinear basis elements at 1 cell-centered
128128
!! quadrature point per cell [L-1 ~> m-1].
129+
real, pointer, dimension(:,:,:) :: Jac => NULL() !< Jacobian determinant |J_q| = a_q*d_q of the element
130+
!! mapping at each of the 4 Gaussian quadrature points [L2 ~> m2].
131+
!! Equal to G%areaT only for rectangular elements; differs when
132+
!! opposite cell edges have unequal lengths (non-rectangular quads).
129133
real, pointer, dimension(:,:,:,:,:,:) :: Phisub => NULL() !< Quadrature structure weights at subgridscale
130134
!! locations for finite element calculations [nondim]
131135
integer :: OD_rt_counter = 0 !< A counter of the number of contributions to OD_rt.
@@ -672,8 +676,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
672676
endif
673677

674678
allocate(CS%Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0)
679+
allocate(CS%Jac(1:4,isd:ied,jsd:jed), source=0.0)
675680
do j=G%jsd,G%jed ; do i=G%isd,G%ied
676-
call bilinear_shape_fn_grid(G, i, j, CS%Phi(:,:,i,j))
681+
call bilinear_shape_fn_grid(G, i, j, CS%Phi(:,:,i,j), CS%Jac(:,i,j))
677682
enddo; enddo
678683

679684
if (CS%GL_regularize) then
@@ -2679,6 +2684,7 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask,
26792684

26802685
real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1]
26812686
real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1]
2687+
real :: jac_wt ! Per-quadrature-point metric correction |J_q|/areaT [nondim]
26822688
integer :: iq, jq, iphi, jphi, i, j, ilq, jlq, Itgt, Jtgt, qp, qpv
26832689
logical :: visc_qp4
26842690
real, dimension(2) :: xquad ! Nondimensional quadrature ratios [nondim]
@@ -2738,22 +2744,25 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask,
27382744
(v_shlf(I-1,J) * Phi(6,qp,i,j)))
27392745

27402746
if (visc_qp4) qpv = qp !current quad point for viscosity
2747+
! Ratio |J_q|/areaT corrects the uniform-area weight baked into ice_visc for
2748+
! non-rectangular elements where opposite cell edges have unequal lengths.
2749+
jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j)
27412750

27422751
do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi
2743-
if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * &
2752+
if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = jac_wt * ice_visc(i,j,qpv) * &
27442753
(((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + &
27452754
((uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)))
2746-
if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * &
2755+
if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = jac_wt * ice_visc(i,j,qpv) * &
27472756
(((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + &
27482757
((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)))
27492758

27502759
if (float_cond(i,j) == 0) then
27512760
ilq = 1 ; if (iq == iphi) ilq = 2
27522761
jlq = 1 ; if (jq == jphi) jlq = 2
27532762
if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + &
2754-
((basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)))
2763+
(jac_wt * (basal_trac(i,j) * uq)) * (xquad(ilq) * xquad(jlq))
27552764
if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + &
2756-
((basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)))
2765+
(jac_wt * (basal_trac(i,j) * vq)) * (xquad(ilq) * xquad(jlq))
27572766
endif
27582767
enddo ; enddo
27592768
enddo ; enddo
@@ -2942,6 +2951,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac,
29422951

29432952
real :: ux, uy, vx, vy ! Interpolated weight gradients [L-1 ~> m-1]
29442953
real :: uq, vq
2954+
real :: jac_wt ! Per-quadrature-point metric correction |J_q|/areaT [nondim]
29452955
real, dimension(2) :: xquad
29462956
real, dimension(2,2) :: Hcell, sub_ground
29472957
real, dimension(2,2,4) :: u_diag_qp, v_diag_qp
@@ -2974,6 +2984,9 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac,
29742984

29752985
qp = 2*(jq-1)+iq !current quad point
29762986
if (visc_qp4) qpv = qp !current quad point for viscosity
2987+
! Ratio |J_q|/areaT corrects the uniform-area weight baked into ice_visc for
2988+
! non-rectangular elements where opposite cell edges have unequal lengths.
2989+
jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j)
29772990

29782991
do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi
29792992

@@ -2987,14 +3000,14 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac,
29873000
vx = 0.
29883001
vy = 0.
29893002

2990-
u_diag_qp(iphi,jphi,qp) = &
3003+
u_diag_qp(iphi,jphi,qp) = jac_wt * &
29913004
ice_visc(i,j,qpv) * (((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + &
29923005
((uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)))
29933006

29943007
if (float_cond(i,j) == 0) then
29953008
uq = xquad(ilq) * xquad(jlq)
29963009
u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + &
2997-
(basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq))
3010+
jac_wt * (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq))
29983011
endif
29993012
endif
30003013

@@ -3005,14 +3018,14 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac,
30053018
ux = 0.
30063019
uy = 0.
30073020

3008-
v_diag_qp(iphi,jphi,qp) = &
3021+
v_diag_qp(iphi,jphi,qp) = jac_wt * &
30093022
ice_visc(i,j,qpv) * (((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + &
30103023
((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)))
30113024

30123025
if (float_cond(i,j) == 0) then
30133026
vq = xquad(ilq) * xquad(jlq)
30143027
v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + &
3015-
(basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq))
3028+
jac_wt * (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq))
30163029
endif
30173030
endif
30183031
enddo ; enddo
@@ -3625,12 +3638,14 @@ end subroutine bilinear_shape_functions
36253638
!> This subroutine calculates the gradients of bilinear basis elements that are centered at the
36263639
!! vertices of the cell using a locally orthogoal MOM6 grid. Values are calculated at
36273640
!! points of gaussian quadrature.
3628-
subroutine bilinear_shape_fn_grid(G, i, j, Phi)
3629-
type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
3641+
subroutine bilinear_shape_fn_grid(G, i, j, Phi, Jac)
3642+
type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
36303643
integer, intent(in) :: i !< The i-index in the grid to work on.
36313644
integer, intent(in) :: j !< The j-index in the grid to work on.
36323645
real, dimension(8,4), intent(inout) :: Phi !< The gradients of bilinear basis elements at Gaussian
36333646
!! quadrature points surrounding the cell vertices [L-1 ~> m-1].
3647+
real, dimension(4), optional, intent(out) :: Jac !< Jacobian determinant |J_q| = a_q*d_q at each
3648+
!! Gaussian quadrature point [L2 ~> m2].
36343649

36353650
! This subroutine calculates the gradients of bilinear basis elements that
36363651
! that are centered at the vertices of the cell. The values are calculated at
@@ -3684,6 +3699,7 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi)
36843699
Phi(2*node,qpoint) = ( (a * (2 * ynode - 3)) * xexp ) / (a*d)
36853700

36863701
enddo
3702+
if (present(Jac)) Jac(qpoint) = a * d
36873703
enddo
36883704

36893705
end subroutine bilinear_shape_fn_grid
@@ -4003,6 +4019,7 @@ subroutine ice_shelf_dyn_end(CS)
40034019
deallocate(CS%OD_rt, CS%OD_av)
40044020
deallocate(CS%t_bdry_val, CS%bed_elev)
40054021
deallocate(CS%ground_frac, CS%ground_frac_rt)
4022+
if (associated(CS%Jac)) deallocate(CS%Jac)
40064023

40074024
deallocate(CS)
40084025

0 commit comments

Comments
 (0)
⚡