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

Add free-slip, partial-slip and no-slip momentum BCs #49

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Add free-slip, partial-slip and no-slip momentum BCs #49

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
74540c0
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 15, 2023
bbd2b9c
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 15, 2023
7947475
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 15, 2023
e9d759a
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 15, 2023
6be1fe7
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 16, 2023
a6f60e1
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 16, 2023
548bdbf
Add a free-slip/partial-slip/no-slip lateral boundary condition
dengwirda Apr 18, 2023
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
6 changes: 6 additions & 0 deletions components/mpas-ocean/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,6 +202,12 @@
possible_values="'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit'"
/>
</nml_record>
<nml_record name="lateral_walls" mode="forward">
<nml_option name="config_wall_slip_factor" type="real" default_value="0.0"
description="Lateral wall slip boundary condition: no-slip=0.0, free-slip=1.0, with partial-slip lying in-between. This factor multiplies relative vorticity computed at vertices positioned on lateral boundaries."
possible_values="Any real number between 0.0 and 1.0"
/>
</nml_record>
<nml_record name="hmix" mode="forward">
<nml_option name="config_hmix_scaleWithMesh" type="logical" default_value=".false."
description="If false, del2 and del4 coefficients are constant throughout the mesh (equivalent to setting $\rho_m=1$ throughout the mesh). If true, these coefficients scale as mesh density to the -3/4 power."
Expand Down
60 changes: 38 additions & 22 deletions components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -567,7 +567,7 @@ end subroutine ocn_diagnostic_solve!}}}
! routine ocn_relativeVorticity_circulation
!
!> \brief Computes relative vorticity and circulation
!> \author Mark Petersen, Doug Jacobsen, Todd Ringler
!> \author Mark Petersen, Doug Jacobsen, Todd Ringler, Darren Engwirda
!> \date November 2013
!> \details
!> Computes relative vorticity and circulation
Expand Down Expand Up @@ -607,7 +607,7 @@ subroutine ocn_relativeVorticity_circulation(relativeVorticity, circulation, nor

integer :: iVertex, iEdge, i, k

real (kind=RKIND) :: invAreaTri1, r_tmp
real (kind=RKIND) :: invAreaTri1, r_tmp, wall_slip_factor

err = 0

Expand All @@ -627,24 +627,36 @@ subroutine ocn_relativeVorticity_circulation(relativeVorticity, circulation, nor
!$omp end do
#endif

! no-slip = 0.0
! free-slip = 1.0
! partial-slip in-between
wall_slip_factor = config_wall_slip_factor

#ifdef MPAS_OPENACC
!$acc parallel loop &
!$acc present(circulation, relativeVorticity, areaTriangle, edgesOnVertex, &
!$acc maxLevelVertexBot, dcEdge, normalVelocity, edgeSignOnVertex, &
!$acc minLevelVertexTop) &
!$acc private(invAreaTri1, iVertex, iEdge, k, r_tmp)
!$acc present(circulation, relativeVorticity, edgesOnVertex, &
!$acc maxLevelVertexBot, dcEdge, normalVelocity, &
!$acc minLevelVertexTop, areaTriangle, edgeSignOnVertex) &
!$acc private(invAreaTri1, i, iEdge, k, r_tmp, wall_slip_factor)
#else
!$omp do schedule(runtime) private(invAreaTri1, i, iEdge, k, r_tmp)
!$omp do schedule(runtime) &
!$omp private(invAreaTri1, i, iEdge, k, r_tmp, wall_slip_factor)
Comment on lines -632 to +643
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is boundaryVertex also needed in the present() arguments?

#endif
do iVertex = 1, nVerticesAll
invAreaTri1 = 1.0_RKIND / areaTriangle(iVertex)
do i = 1, vertexDegree
iEdge = edgesOnVertex(i, iVertex)
do k = minLevelVertexTop(iVertex), maxLevelVertexBot(iVertex)
r_tmp = dcEdge(iEdge) * normalVelocity(k, iEdge)
circulation(k, iVertex) = circulation(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp
relativeVorticity(k, iVertex) = relativeVorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp * invAreaTri1
do k = minLevelVertexTop(iVertex), maxLevelVertexBot(iVertex)
do i = 1, vertexDegree
iEdge = edgesOnVertex(i, iVertex)
r_tmp = edgeSignOnVertex(i, iVertex) * &
dcEdge(iEdge) * normalVelocity(k, iEdge)
circulation(k, iVertex) = circulation(k, iVertex) + r_tmp
end do
! no-slip: curl(u) unchanged
! free-slip: curl(u) = 0.0
! partial-slip: curl(u) reduced
circulation(k, iVertex) = circulation(k, iVertex) * &
(1.0_RKIND - wall_slip_factor * boundaryVertex(k, iVertex))
relativeVorticity(k, iVertex) = circulation(k, iVertex) * invAreaTri1
Copy link

Choose a reason for hiding this comment

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

The changes in this loop are what is needed to include wall_slip_factor. For this section, it is true that wall_slip_factor = 0.0 should only change the solution due to order of operations. For a second PR for just the addition of wall_slip_factor and these lines, it will be easier to show that only minor changes occurred.

Copy link
Collaborator

Choose a reason for hiding this comment

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

@dengwirda, is there a reason to exchange loop order here? It seems like the order of operations issue could be avoided by applying the (1.0_RKIND - wall_slip_factor * boundaryVertex(k, iVertex)) factor in a separate loop k loop following the original do i = 1, vertexDegree loop.

end do
end do
#ifndef MPAS_OPENACC
Expand Down Expand Up @@ -1008,11 +1020,11 @@ subroutine ocn_diagnostic_solve_vorticity(dt, &
cell1, cell2 ! neighbor cell addresses

real(kind=RKIND) :: &
invAreaTri1, &! 1/triangle area
invAreaCell1, &! 1/cell area
invLength, &! 1/length
layerThicknessVertex, &! layer thickness at vertex
apvm_scale_factor ! scale factor for APVM form
apvm_scale_factor, &! scale factor for APVM form
areaDual ! unmaksed dual-cell overlap

! Local arrays
! normalizedPlanetaryVorticityVertex: earth's rotational rate
Expand Down Expand Up @@ -1052,25 +1064,29 @@ subroutine ocn_diagnostic_solve_vorticity(dt, &
!$acc present(normalizedRelativeVorticityVertex, &
!$acc normalizedPlanetaryVorticityVertex, &
!$acc relativeVorticity, fVertex, layerThickness, &
!$acc areaTriangle, kiteAreasOnVertex, cellsOnVertex, &
!$acc maxLevelVertexBot) &
!$acc private(i, k, iCell, invAreaTri1, layerThicknessVertex)
!$acc kiteAreasOnVertex, cellsOnVertex, &
!$acc minLevelVertexTop, maxLevelVertexBot) &
!$acc private(i, k, iCell, areaDual, layerThicknessVertex)
#else
!$omp parallel
!$omp do schedule(runtime) &
!$omp private(i, k, iCell, invAreaTri1, layerThicknessVertex)
!$omp private(i, k, iCell, areaDual, layerThicknessVertex)
#endif
do iVertex = 1, nVertices
invAreaTri1 = 1.0_RKIND / areaTriangle(iVertex)
do k = 1, maxLevelVertexBot(iVertex)
normalizedRelativeVorticityVertex(:,iVertex) = 0.0_RKIND
normalizedPlanetaryVorticityVertex(:,iVertex) = 0.0_RKIND
do k = minLevelVertexTop(iVertex), maxLevelVertexBot(iVertex)
layerThicknessVertex = 0.0_RKIND
areaDual = 0.0_RKIND ! only the unmasked kites
do i = 1, vertexDegree
iCell = cellsOnVertex(i,iVertex)
layerThicknessVertex = layerThicknessVertex &
+ layerThickness(k,iCell) &
* kiteAreasOnVertex(i,iVertex)
areaDual = areaDual + cellMask(k,iCell) &
* kiteAreasOnVertex(i,iVertex)
end do
layerThicknessVertex = layerThicknessVertex*invAreaTri1
layerThicknessVertex = layerThicknessVertex/areaDual
Copy link

Choose a reason for hiding this comment

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

@xylar I see what you are saying. This is the change where invAreaTri1 includes all three kites, while areaDual only includes the kites from active cells. This would be answer changing, and the changes in the loop above could go into a single PR.

if (layerThicknessVertex == 0) cycle

normalizedRelativeVorticityVertex(k,iVertex) = &
Expand Down