From 77e22f0edace70ab58064108f8a81da950c11e58 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 22 Jan 2021 12:14:01 -0800 Subject: [PATCH 1/4] Do not zero N under floating ice in Interface This commit removes code that zeros effective pressure beneath floating ice in the Interface to Albany. This is unneeded because Albany zeros beta at quadrature points, and zeroing it here at FEM nodes interferes with that. At best it's redundant, and at worst (which is if N were to be assigned nonzero values under floating ice), it results in significant errors in the solve. This is also redundant because MPAS hydro model and ocean_connection_N parameterization will already zero N under floating ice. --- src/core_landice/mode_forward/Interface_velocity_solver.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index 3971076ba9..d72df8476c 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -1252,8 +1252,6 @@ void importFields(std::vector >& marineBdyExtensionMap, dou } bool floating = rho_ice * thicknessData[iv] + rho_ocean * bedTopographyData[iv] < 0; - if (floating && (effecPress_F != 0)) - effecPressData[iv] = 0.0; } } From cdfee640432e4fa896fff881db5bb33b9f8c2350 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 27 Jan 2021 08:03:18 -0800 Subject: [PATCH 2/4] Add option to extend effective pressure for the velocity solver This is needed because the friction laws that use effective pressure need N on the first floating cell to allow for more accurate evaluation of friction at quadrature points. Also, this commit moves the limiting of the effective pressure from the velocity_external routine to the calc_beta routine. Finally, all these manipulations to N are done on the field effectivePressureLimited to avoid messing up effectivePressure which is used by the subglacial hydrology model (if on). --- src/core_landice/Registry.xml | 4 ++ .../mode_forward/mpas_li_velocity.F | 44 +++++++++++++++++++ .../mode_forward/mpas_li_velocity_external.F | 9 ---- 3 files changed, 48 insertions(+), 9 deletions(-) diff --git a/src/core_landice/Registry.xml b/src/core_landice/Registry.xml index af70c9f57d..8ba24d5831 100644 --- a/src/core_landice/Registry.xml +++ b/src/core_landice/Registry.xml @@ -100,6 +100,10 @@ description="If true, nonconvergence of the velocity solver results in a fatal error. If false, it is a warning." possible_values=".true. or .false." /> + 0) then + effectivePressureLimited(iCell) = NAccum / real(nNValues, kind=RKIND) ! take average of neighboring effectivePressure values + ! < TODO: could do this in log space + endif + endif ! if iCell is floating + enddo ! cell loop + endif + + ! Apply limits to N + do iCell = 1, nCells + effectivePressureLimited(iCell) = max(0.0_RKIND, effectivePressureLimited(iCell)) ! Never pass a negative N + effectivePressureLimited(iCell) = min(config_effective_pressure_max, effectivePressureLimited(iCell)) + enddo + + !-------------------------------------------------------------------- end subroutine calculate_beta diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F index 10d3b38914..519dbeb66c 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity_external.F +++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F @@ -449,7 +449,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro normalVelocity, uReconstructX, uReconstructY, uReconstructZ real (kind=RKIND), dimension(:,:), pointer :: temperature real (kind=RKIND), dimension(:), pointer :: stiffnessFactor - real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited real (kind=RKIND), dimension(:), pointer :: muFriction real (kind=RKIND), pointer :: deltat @@ -460,7 +459,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro logical, pointer :: config_output_external_velocity_solver_data logical, pointer :: config_nonconvergence_error real (kind=RKIND), pointer :: config_ice_density - real (kind=RKIND), pointer :: config_effective_pressure_max integer, pointer :: anyDynamicVertexMaskChanged integer, pointer :: dirichletMaskChanged integer, pointer :: nEdges @@ -486,7 +484,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro config_output_external_velocity_solver_data) call mpas_pool_get_config(liConfigs, 'config_nonconvergence_error', config_nonconvergence_error) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) - call mpas_pool_get_config(liConfigs, 'config_effective_pressure_max', config_effective_pressure_max) ! Mesh variables call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) @@ -524,7 +521,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) ! Hydro variables - call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure) call mpas_pool_get_array(velocityPool, 'effectivePressureLimited', effectivePressureLimited) #if defined(USE_EXTERNAL_L1L2) || defined(USE_EXTERNAL_FIRSTORDER) || defined(USE_EXTERNAL_STOKES) @@ -532,11 +528,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call interface_redirect_stdout(timestepNumber) #endif - ! Create version of effectivePressure field with limits applied - do iCell = 1, nCells - effectivePressureLimited(iCell) = max(0.0_RKIND, effectivePressure(iCell)) ! Never pass a negative N - effectivePressureLimited(iCell) = min(config_effective_pressure_max, effectivePressureLimited(iCell)) - enddo ! ================================================================== ! External dycore calls to be made only when vertex mask changes From c4805269e208952dd0b487e60d525b9cdbaf0fec Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 29 Jan 2021 12:10:22 -0800 Subject: [PATCH 3/4] Change from mean to min extrap of N --- src/core_landice/mode_forward/mpas_li_velocity.F | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/core_landice/mode_forward/mpas_li_velocity.F b/src/core_landice/mode_forward/mpas_li_velocity.F index eeec864edc..d87163d35a 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity.F +++ b/src/core_landice/mode_forward/mpas_li_velocity.F @@ -819,6 +819,7 @@ subroutine calculate_beta(block, err) real (kind=RKIND) :: betaAccum integer :: nBetaValues real (kind=RKIND) :: NAccum + real (kind=RKIND) :: Nmin integer :: nNValues integer :: iCell, iCell2 integer :: neighbor @@ -868,16 +869,19 @@ subroutine calculate_beta(block, err) if ( .not. li_mask_is_grounded_ice(cellMask(iCell)) ) then nNValues = 0 NAccum = 0.0_RKIND + Nmin = 1.0e20 do iCell2 = 1, nEdgesOnCell(iCell) neighbor = cellsOnCell(iCell2, iCell) if (li_mask_is_grounded_ice(cellMask(neighbor))) then NAccum = NAccum + effectivePressureLimited(neighbor) nNValues = nNValues + 1 + Nmin = min(Nmin, effectivePressureLimited(neighbor)) endif enddo if (nNValues > 0) then - effectivePressureLimited(iCell) = NAccum / real(nNValues, kind=RKIND) ! take average of neighboring effectivePressure values +! effectivePressureLimited(iCell) = NAccum / real(nNValues, kind=RKIND) ! take average of neighboring effectivePressure values ! < TODO: could do this in log space + effectivePressureLimited(iCell) = Nmin endif endif ! if iCell is floating enddo ! cell loop From 8650558db7b4509201a843772773f37059f475f8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 4 Mar 2021 07:30:19 -0800 Subject: [PATCH 4/4] HACK: Make N=mu*N to allow N-based runs This hack sets N=mu*N before passing N to Albany. This is hack to get Albany to only evaluate mu*N at cell centers rather than at all quadrature points. This hack is used for doing a forward run where mu is calculated from an MPAS N field. In this configuration, the albany input file should mu be spatially uniform with a value of 1.0. This commit should be removed or formalized before this branch is merged!!! --- Makefile | 52 +++++++++++++++++++ .../mode_forward/mpas_li_velocity.F | 3 ++ 2 files changed, 55 insertions(+) diff --git a/Makefile b/Makefile index 296013999a..ddc04acea5 100644 --- a/Makefile +++ b/Makefile @@ -454,6 +454,58 @@ intel-nersc: "OPENMP = $(OPENMP)" \ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" ) +intel-cori-knl: + ( $(MAKE) all \ + "FC_PARALLEL = ftn" \ + "CC_PARALLEL = cc" \ + "CXX_PARALLEL = CC" \ + "FC_SERIAL = ftn" \ + "CC_SERIAL = cc" \ + "CXX_SERIAL = CC" \ + "FFLAGS_PROMOTION = -real-size 64" \ + "FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \ + "CFLAGS_OPT = -O3" \ + "CXXFLAGS_OPT = -O3, -xMIC-AVX512" \ + "LDFLAGS_OPT = -O3" \ + "FFLAGS_OMP = -qopenmp" \ + "CFLAGS_OMP = -qopenmp" \ + "FFLAGS_DEBUG = -real-size 64 -g -convert big_endian -free -CU -CB -check all -gen-interfaces -warn interfaces -traceback" \ + "CFLAGS_DEBUG = -g -traceback" \ + "CXXFLAGS_DEBUG = -g -traceback, -xMIC-AVX512" \ + "LDFLAGS_DEBUG = -g -traceback" \ + "BUILD_TARGET = $(@)" \ + "CORE = $(CORE)" \ + "DEBUG = $(DEBUG)" \ + "USE_PAPI = $(USE_PAPI)" \ + "OPENMP = $(OPENMP)" \ + "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" ) + +intel-cori-haswell: + ( $(MAKE) all \ + "FC_PARALLEL = ftn" \ + "CC_PARALLEL = cc" \ + "CXX_PARALLEL = CC" \ + "FC_SERIAL = ftn" \ + "CC_SERIAL = cc" \ + "CXX_SERIAL = CC" \ + "FFLAGS_PROMOTION = -real-size 64" \ + "FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \ + "CFLAGS_OPT = -O3" \ + "CXXFLAGS_OPT = -O3, -xCORE-AVX2" \ + "LDFLAGS_OPT = -O3" \ + "FFLAGS_OMP = -qopenmp" \ + "CFLAGS_OMP = -qopenmp" \ + "FFLAGS_DEBUG = -real-size 64 -g -convert big_endian -free -CU -CB -check all -gen-interfaces -warn interfaces -traceback" \ + "CFLAGS_DEBUG = -g -traceback" \ + "CXXFLAGS_DEBUG = -g -traceback, -xCORE-AVX2" \ + "LDFLAGS_DEBUG = -g -traceback" \ + "BUILD_TARGET = $(@)" \ + "CORE = $(CORE)" \ + "DEBUG = $(DEBUG)" \ + "USE_PAPI = $(USE_PAPI)" \ + "OPENMP = $(OPENMP)" \ + "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" ) + bluegene: ( $(MAKE) all \ "FC_PARALLEL = mpixlf95_r" \ diff --git a/src/core_landice/mode_forward/mpas_li_velocity.F b/src/core_landice/mode_forward/mpas_li_velocity.F index d87163d35a..76a7f8b157 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity.F +++ b/src/core_landice/mode_forward/mpas_li_velocity.F @@ -810,6 +810,7 @@ subroutine calculate_beta(block, err) real (kind=RKIND), dimension(:), pointer :: betaSolve, beta real (kind=RKIND), dimension(:), pointer :: effectivePressure, effectivePressureLimited real (kind=RKIND), dimension(:), pointer :: basalTemperature, basalPmpTemperature + real (kind=RKIND), dimension(:), pointer :: muFriction integer, dimension(:), pointer :: cellMask logical, pointer :: config_use_glp logical, pointer :: config_beta_thawed_only @@ -842,6 +843,7 @@ subroutine calculate_beta(block, err) call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve) call mpas_pool_get_array(velocityPool, 'beta', beta) + call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure) call mpas_pool_get_array(velocityPool, 'effectivePressureLimited', effectivePressureLimited) @@ -892,6 +894,7 @@ subroutine calculate_beta(block, err) effectivePressureLimited(iCell) = max(0.0_RKIND, effectivePressureLimited(iCell)) ! Never pass a negative N effectivePressureLimited(iCell) = min(config_effective_pressure_max, effectivePressureLimited(iCell)) enddo + effectivePressureLimited(:) = effectivePressureLimited(:) * muFriction(:) !--------------------------------------------------------------------