Skip to content

Commit

Permalink
Merge branch 'akturner/mpas-seaice/refactor-velocity-init' (PR #5290)
Browse files Browse the repository at this point in the history 
Refactor of MPAS-Seaice velocity solver

Refactor of MPAS sea ice velocity solver to allow easier testing of new
options and methods, including:
* Move domain reference to within subroutines avoiding messy passing of
  variables
* Move wachspress basis routines to separate module
* Move triangular quadrature rules to new module
* Move some mesh related routines to the mesh module
* Rename block variable to blockPtr
* Rename all pool types to end in Pool

Fixes #5289

[BFB]
  • Loading branch information
@jonbob
jonbob committed Dec 15, 2022
2 parents 83d265b + 99715da commit 4b872d5
Show file tree
Hide file tree
Showing 12 changed files with 3,535 additions and 3,556 deletions.
2 changes: 2 additions & 0 deletions components/mpas-seaice/src/seaice.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,8 @@ list(APPEND RAW_SOURCES
core_seaice/shared/mpas_seaice_velocity_solver_pwl.F
core_seaice/shared/mpas_seaice_velocity_solver_variational_shared.F
core_seaice/shared/mpas_seaice_velocity_solver_constitutive_relation.F
core_seaice/shared/mpas_seaice_triangle_quadrature.F
core_seaice/shared/mpas_seaice_wachspress_basis.F
core_seaice/shared/mpas_seaice_forcing.F
core_seaice/shared/mpas_seaice_initialize.F
core_seaice/shared/mpas_seaice_testing.F
Expand Down
10 changes: 8 additions & 2 deletions components/mpas-seaice/src/shared/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ OBJS = mpas_seaice_time_integration.o \
mpas_seaice_velocity_solver_pwl.o \
mpas_seaice_velocity_solver_variational_shared.o \
mpas_seaice_velocity_solver_constitutive_relation.o \
mpas_seaice_wachspress_basis.o \
mpas_seaice_triangle_quadrature.o \
mpas_seaice_forcing.o \
mpas_seaice_initialize.o \
mpas_seaice_testing.o \
Expand Down Expand Up @@ -48,15 +50,19 @@ mpas_seaice_velocity_solver_constitutive_relation.o: mpas_seaice_constants.o mpa

mpas_seaice_forcing.o: mpas_seaice_constants.o mpas_seaice_mesh.o

mpas_seaice_wachspress_basis.o: mpas_seaice_mesh.o

mpas_seaice_triangle_quadrature.o:

mpas_seaice_velocity_solver_weak.o: mpas_seaice_constants.o mpas_seaice_testing.o mpas_seaice_velocity_solver_constitutive_relation.o mpas_seaice_mesh_pool.o

mpas_seaice_velocity_solver_variational_shared.o: mpas_seaice_constants.o

mpas_seaice_velocity_solver_wachspress.o: mpas_seaice_constants.o mpas_seaice_numerics.o mpas_seaice_mesh.o mpas_seaice_testing.o mpas_seaice_velocity_solver_variational_shared.o
mpas_seaice_velocity_solver_wachspress.o: mpas_seaice_constants.o mpas_seaice_numerics.o mpas_seaice_mesh.o mpas_seaice_testing.o mpas_seaice_velocity_solver_variational_shared.o mpas_seaice_wachspress_basis.o mpas_seaice_triangle_quadrature.o

mpas_seaice_velocity_solver_pwl.o: mpas_seaice_constants.o mpas_seaice_numerics.o mpas_seaice_mesh.o mpas_seaice_testing.o mpas_seaice_velocity_solver_variational_shared.o

mpas_seaice_velocity_solver_variational.o: mpas_seaice_constants.o mpas_seaice_velocity_solver_constitutive_relation.o mpas_seaice_velocity_solver_wachspress.o mpas_seaice_velocity_solver_pwl.o mpas_seaice_mesh_pool.o
mpas_seaice_velocity_solver_variational.o: mpas_seaice_constants.o mpas_seaice_velocity_solver_constitutive_relation.o mpas_seaice_velocity_solver_wachspress.o mpas_seaice_velocity_solver_pwl.o mpas_seaice_mesh_pool.o mpas_seaice_mesh.o

mpas_seaice_velocity_solver.o: mpas_seaice_constants.o mpas_seaice_mesh.o mpas_seaice_testing.o mpas_seaice_velocity_solver_weak.o mpas_seaice_velocity_solver_constitutive_relation.o mpas_seaice_velocity_solver_variational.o mpas_seaice_diagnostics.o mpas_seaice_mesh_pool.o mpas_seaice_special_boundaries.o

Expand Down
269 changes: 238 additions & 31 deletions components/mpas-seaice/src/shared/mpas_seaice_mesh.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,8 @@ module seaice_mesh

public :: &
seaice_init_mesh, &
seaice_cell_vertices_at_vertex, &
seaice_wrapped_index, &
seaice_calc_local_coords, &
seaice_normal_vectors, &
seaice_normal_vectors_polygon, &
seaice_dot_product_3space, &
Expand Down Expand Up @@ -614,12 +615,12 @@ subroutine interior_edges(&
end subroutine interior_edges

!-----------------------------------------------------------------------
! mesh searches
! misc
!-----------------------------------------------------------------------

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! seaice_cell_vertices_at_vertex
! seaice_wrapped_index
!
!> \brief
!> \author Adrian K. Turner, LANL
Expand All @@ -629,60 +630,266 @@ end subroutine interior_edges
!
!-----------------------------------------------------------------------

subroutine seaice_cell_vertices_at_vertex(&
cellVerticesAtVertex, &
nVertices, &
vertexDegree, &
function seaice_wrapped_index(&
input, &
nelements) &
result(output)!{{{

integer, intent(in) :: &
input, & !< Input:
nelements !< Input:

integer :: output

output = modulo(input - 1, nelements) + 1

end function seaice_wrapped_index!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! seaice_calc_local_coords
!
!> \brief
!> \author Adrian K. Turner, LANL
!> \date 22 October 2014
!> \details
!>
!
!-----------------------------------------------------------------------

subroutine seaice_calc_local_coords(&
xLocal, &
yLocal, &
nCells, &
nEdgesOnCell, &
verticesOnCell, &
xVertex, &
yVertex, &
zVertex, &
xCell, &
yCell, &
zCell, &
rotateCartesianGrid, &
onASphere)!{{{

real(kind=RKIND), dimension(:,:), intent(out) :: &
xLocal, & !< Output:
yLocal !< Output:

integer, intent(in) :: &
nCells !< Input:

integer, dimension(:), intent(in) :: &
nEdgesOnCell !< Input:

integer, dimension(:,:), intent(in) :: &
verticesOnCell !< Input:

real(kind=RKIND), dimension(:), intent(in) :: &
xVertex, & !< Input:
yVertex, & !< Input:
zVertex, & !< Input:
xCell, & !< Input:
yCell, & !< Input:
zCell !< Input:

logical, intent(in) :: &
rotateCartesianGrid, & !< Input:
onASphere !< Input:

if (onASphere) then
call calc_local_coords_spherical(&
xLocal, &
yLocal, &
nCells, &
nEdgesOnCell, &
verticesOnCell, &
xVertex, &
yVertex, &
zVertex, &
xCell, &
yCell, &
zCell, &
rotateCartesianGrid)
else
call calc_local_coords_planar(&
xLocal, &
yLocal, &
nCells, &
nEdgesOnCell, &
verticesOnCell, &
xVertex, &
yVertex, &
zVertex, &
xCell, &
yCell, &
zCell)
endif

end subroutine seaice_calc_local_coords!}}}

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! calc_local_coords_planar
!
!> \brief
!> \author Adrian K. Turner, LANL
!> \date 2013-2014
!> \details
!>
!
!-----------------------------------------------------------------------

subroutine calc_local_coords_planar(&
xLocal, &
yLocal, &
nCells, &
nEdgesOnCell, &
verticesOnCell, &
cellsOnVertex)!{{{
xVertex, &
yVertex, &
zVertex, &
xCell, &
yCell, &
zCell)!{{{

integer, dimension(:,:), intent(out) :: &
cellVerticesAtVertex !< Output:
real(kind=RKIND), dimension(:,:), intent(out) :: &
xLocal, & !< Output:
yLocal !< Output:

integer, intent(in) :: &
nVertices, & !< Input:
vertexDegree !< Input:
nCells !< Input:

integer, dimension(:), intent(in) :: &
nEdgesOnCell !< Input:

integer, dimension(:,:), intent(in) :: &
cellsOnVertex, & !< Input:
verticesOnCell !< Input:
verticesOnCell !< Input:

real(kind=RKIND), dimension(:), intent(in) :: &
xVertex, & !< Input:
yVertex, & !< Input:
zVertex, & !< Input:
xCell, & !< Input:
yCell, & !< Input:
zCell !< Input:

integer :: &
iVertex, &
iVertexDegree, &
iCell, &
iVertexOnCell, &
jVertex
iVertex, &
iVertexOnCell

do iVertex = 1, nVertices
do iCell = 1, nCells

do iVertexDegree = 1, vertexDegree
do iVertexOnCell = 1, nEdgesOnCell(iCell)

cellVerticesAtVertex(iVertexDegree,iVertex) = 0
iVertex = verticesOnCell(iVertexOnCell, iCell)

iCell = cellsOnVertex(iVertexDegree, iVertex)
xLocal(iVertexOnCell,iCell) = xVertex(iVertex) - xCell(iCell)
yLocal(iVertexOnCell,iCell) = yVertex(iVertex) - yCell(iCell)

do iVertexOnCell = 1, nEdgesOnCell(iCell)
enddo ! iVertexOnCell

jVertex = verticesOnCell(iVertexOnCell,iCell)
enddo ! iCell

if (iVertex == jVertex) then
end subroutine calc_local_coords_planar!}}}

cellVerticesAtVertex(iVertexDegree,iVertex) = iVertexOnCell
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! calc_local_coords_spherical
!
!> \brief
!> \author Adrian K. Turner, LANL
!> \date 22 October 2014
!> \details
!>
!
!-----------------------------------------------------------------------

endif
subroutine calc_local_coords_spherical(&
xLocal, &
yLocal, &
nCells, &
nEdgesOnCell, &
verticesOnCell, &
xVertex, &
yVertex, &
zVertex, &
xCell, &
yCell, &
zCell, &
rotateCartesianGrid)!{{{

real(kind=RKIND), dimension(:,:), intent(out) :: &
xLocal, & !< Output:
yLocal !< Output:

enddo ! iVertexOnCell
integer, intent(in) :: &
nCells !< Input:

enddo ! iVertexDegree
integer, dimension(:), intent(in) :: &
nEdgesOnCell !< Input:

enddo ! iVertex
integer, dimension(:,:), intent(in) :: &
verticesOnCell !< Input:

real(kind=RKIND), dimension(:), intent(in) :: &
xVertex, & !< Input:
yVertex, & !< Input:
zVertex, & !< Input:
xCell, & !< Input:
yCell, & !< Input:
zCell !< Input:

logical, intent(in) :: &
rotateCartesianGrid !< Input:

real(kind=RKIND), dimension(3) :: &
normalVector3D

real(kind=RKIND), dimension(2) :: &
normalVector2D

integer :: &
iCell, &
iVertex, &
iVertexOnCell

real(kind=RKIND) :: &
xCellRotated, &
yCellRotated, &
zCellRotated

do iCell = 1, nCells

call seaice_grid_rotation_forward(&
xCellRotated, yCellRotated, zCellRotated, &
xCell(iCell), yCell(iCell), zCell(iCell), &
rotateCartesianGrid)

do iVertexOnCell = 1, nEdgesOnCell(iCell)

iVertex = verticesOnCell(iVertexOnCell, iCell)

call seaice_grid_rotation_forward(&
normalVector3D(1), normalVector3D(2), normalVector3D(3), &
xVertex(iVertex), yVertex(iVertex), zVertex(iVertex), &
rotateCartesianGrid)

call seaice_project_3D_vector_onto_local_2D(&
normalVector2D, &
normalVector3D, &
xCellRotated, &
yCellRotated, &
zCellRotated)

xLocal(iVertexOnCell,iCell) = normalVector2D(1)
yLocal(iVertexOnCell,iCell) = normalVector2D(2)

enddo ! iVertexOnCell

enddo ! iCell

end subroutine seaice_cell_vertices_at_vertex!}}}
end subroutine calc_local_coords_spherical!}}}

!-----------------------------------------------------------------------
! normal vectors
Expand Down
Loading

0 comments on commit 4b872d5

Please sign in to comment.