From fa6b5ca0bcccdcd613e7c0a6edf23175c3fcb7f9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 28 Oct 2024 19:52:52 -0700 Subject: [PATCH 1/4] Add new spatial bed roughness variable for reg. Coulomb friction This commit adds a new edRoughnessRC variable to MALI on the MPAS side. Corresponding changes are needed on the Albany side. --- components/mpas-albany-landice/src/Registry.xml | 3 +++ .../mode_forward/Interface_velocity_solver.cpp | 17 +++++++++++++---- .../mode_forward/Interface_velocity_solver.hpp | 3 ++- .../mode_forward/mpas_li_velocity_external.F | 6 ++++++ 4 files changed, 24 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 3c4b562c4e3d..6a12649a12cc 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1638,6 +1638,9 @@ is the value of that variable from the *previous* time level! + diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp index 8ecff9880441..357ffec1577f 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp @@ -69,7 +69,7 @@ std::vector indexToTriangleID, std::vector indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge, fVertexToTriangleID, fCellToVertex, iceMarginEdgesLIds, dirichletNodesIDs; std::vector dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells, - elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell; + elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell; std::vector isVertexBoundary, isBoundaryEdge; // only needed for creating ASCII mesh @@ -241,6 +241,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double* const dirichletVelocityXValue, double* const dirichletVelocitYValue, double* u_normal_F, double* bodyForce_F, double* dissipation_heat_F, double* xVelocityOnCell, double* yVelocityOnCell, double const* deltat, @@ -269,7 +270,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower if (!isDomainEmpty) { std::vector > marineBdyExtensionMap; - importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, temperature_F, smb_F, minThickness); + importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRougnessRC_F, temperature_F, smb_F, minThickness); std::vector regulThk(thicknessData); for (int index = 0; index < nVertices; index++) @@ -286,7 +287,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta, regulThk, levelsNormalizedThickness, elevationData, thicknessData, betaData, bedTopographyData, smbData, - stiffnessFactorData, effecPressData, muFrictionData, + stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData temperatureDataOnPrisms, bodyForceOnBasalCell, dissipationHeatOnPrisms, velocityOnVertices, albany_error, dt); *error=albany_error; @@ -1117,6 +1118,7 @@ double signedTriangleAreaOnSphere(const double* x, const double* y, void importFields(std::vector >& marineBdyExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F, double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const * temperature_F, double const * smb_F, double eps) { int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1; @@ -1131,6 +1133,8 @@ void importFields(std::vector >& marineBdyExtensionMap, dou effecPressData.assign(nVertices, 1e10); if (muFriction_F!= 0) muFrictionData.assign(nVertices, 1e10); + if (bedRoughnessRC_F!= 0) + bedRoughnessData.assign(nVertices, 1e10); if(temperature_F != 0) temperatureDataOnPrisms.assign(nLayers * nTriangles, 1e10); if (smb_F != 0) @@ -1153,6 +1157,8 @@ void importFields(std::vector >& marineBdyExtensionMap, dou effecPressData[index] = effecPress_F[iCell] / unit_length; if (muFriction_F != 0) muFrictionData[index] = muFriction_F[iCell]; + if (bedRoughnessRC_F != 0) + bedRoughnessData[index] = bedRoughnessRC_F[iCell]; } int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles; @@ -1630,6 +1636,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep double const* surfaceAirTemperature_F, double const* basalHeatFlux_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const* thickness_F, double const* thicknessUncertainty_F, double const* smb_F, double const* smbUncertainty_F, double const* bmb_F, double const* bmbUncertainty_F, @@ -1715,7 +1722,8 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep std::vector > marineBdyExtensionMap; // local map to be created by importFields importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, - stiffnessFactor_F, effecPress_F, muFriction_F, temperature_F, smb_F, minThickness); + stiffnessFactor_F, effecPress_F, muFriction_F, bedRougnessRC_F, + temperature_F, smb_F, minThickness); import2DFieldsObservations(marineBdyExtensionMap, thicknessUncertainty_F, @@ -1760,6 +1768,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep write_ascii_mesh_field(betaData, "beta"); write_ascii_mesh_field(muFrictionData, "mu"); write_ascii_mesh_field(muLogData, "mu_log"); + write_ascii_mesh_field(bedRoughnessData, "bed_roughness"); write_ascii_mesh_field(stiffnessFactorData, "stiffening_factor"); write_ascii_mesh_field(stiffnessLogData, "stiffening_factor_log"); diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp index fb8a907216ac..7f27e9aab6cc 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp @@ -90,7 +90,7 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F, void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, - double const* effecPress_F, double const* muFriction_F, + double const* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, double* const dirichletVelocityXValue = 0, double* const dirichletVelocitYValue = 0, double* u_normal_F = 0, double* bodyForce_F = 0, double* dissipation_heat_F = 0, double* xVelocityOnCell = 0, double* yVelocityOnCell = 0, double const * deltat = 0, @@ -134,6 +134,7 @@ void write_ascii_mesh(int const* indexToCellID_F, double const* surfaceAirTemperature_F, double const* basalHeatFlux_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const* thickness_F, double const* thicknessUncertainty_F, double const* smb_F, double const* smbUncertainty_F, double const* bmb_F, double const* bmbUncertainty_F, diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F index 077ff4c1773f..b5eb319655cd 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F @@ -452,6 +452,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited real (kind=RKIND), dimension(:), pointer :: muFriction + real (kind=RKIND), dimension(:), pointer :: bedRoughnessRC real (kind=RKIND), pointer :: deltat integer, dimension(:), pointer :: vertexMask, cellMask, edgeMask integer, dimension(:,:), pointer :: dirichletVelocityMask @@ -516,6 +517,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(velocityPool, 'drivingStressVert', drivingStressVert) call mpas_pool_get_array(velocityPool, 'drivingStress', drivingStress) call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) + call mpas_pool_get_array(velocityPool, 'bedRoughnessRC', bedRoughnessRC) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) @@ -598,6 +600,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, & betaSolve, sfcMassBal, temperature, stiffnessFactor, & effectivePressureLimited, muFriction, & + bedRoughnessRC, & uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1 normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values deltat, albanyVelocityError) ! return values @@ -787,6 +790,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) real (kind=RKIND), dimension(:), pointer :: stiffnessFactor real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: muFriction + real (kind=RKIND), dimension(:), pointer :: bedRoughnessRC integer, dimension(:,:), pointer :: dirichletVelocityMask type (mpas_pool_type), pointer :: meshPool, geometryPool, thermalPool, observationsPool, velocityPool, scratchPool, hydroPool real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density @@ -848,6 +852,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! Velocity variables call mpas_pool_get_array(velocityPool, 'beta', beta) call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) + call mpas_pool_get_array(velocityPool, 'bedRoughnessRC', bedRoughnessRC) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) @@ -901,6 +906,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) surfaceAirTemperature, basalHeatFlux, & stiffnessFactor, & effectivePressure, muFriction, & + bedRoughnessRC, & thickness, thicknessUncertainty, & sfcMassBal, sfcMassBalUncertainty, & floatingBasalMassBal, floatingBasalMassBalUncertainty, & From 1c011891dde70dfafe556d0517aa4de69984e9ed Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 1 Nov 2024 15:56:02 -0600 Subject: [PATCH 2/4] Fixing typos and exporting effective pressure field --- .../mpas-albany-landice/src/Registry.xml | 1 + .../Interface_velocity_solver.cpp | 24 +++++++++++++++---- .../Interface_velocity_solver.hpp | 8 ++++--- 3 files changed, 25 insertions(+), 8 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 6a12649a12cc..3764768ba770 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -826,6 +826,7 @@ + diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp index 357ffec1577f..dcc3192f64c0 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp @@ -240,7 +240,7 @@ void velocity_solver_init_fo(double const *levelsRatio_F) { void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, - double const* effecPress_F, double const* muFriction_F, + double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, double* const dirichletVelocityXValue, double* const dirichletVelocitYValue, double* u_normal_F, double* bodyForce_F, double* dissipation_heat_F, @@ -270,7 +270,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower if (!isDomainEmpty) { std::vector > marineBdyExtensionMap; - importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRougnessRC_F, temperature_F, smb_F, minThickness); + importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, temperature_F, smb_F, minThickness); std::vector regulThk(thicknessData); for (int index = 0; index < nVertices; index++) @@ -287,7 +287,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta, regulThk, levelsNormalizedThickness, elevationData, thicknessData, betaData, bedTopographyData, smbData, - stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData + stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, temperatureDataOnPrisms, bodyForceOnBasalCell, dissipationHeatOnPrisms, velocityOnVertices, albany_error, dt); *error=albany_error; @@ -300,6 +300,8 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower } exportBeta(beta_F); + exportEffectivePressure(effecPress_F); + mapVerticesToCells(velocityOnVertices, &velocityOnCells[0], 2, nLayers, Ordering); @@ -1158,7 +1160,7 @@ void importFields(std::vector >& marineBdyExtensionMap, dou if (muFriction_F != 0) muFrictionData[index] = muFriction_F[iCell]; if (bedRoughnessRC_F != 0) - bedRoughnessData[index] = bedRoughnessRC_F[iCell]; + bedRoughnessData[index] = bedRoughnessRC_F[iCell] / unit_length; } int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles; @@ -1367,6 +1369,18 @@ void exportBeta(double * beta_F) { allToAll (beta_F, sendCellsList_F, recvCellsList_F, 1); } +void exportEffectivePressure(double * effecPress_F) { + std::fill(effecPress_F, effecPress_F + nCells_F, 0.); + for (int index = 0; index < nVertices; index++) { + int fCell = vertexToFCell[index]; + effecPress_F[fCell] = effecPressData[index] * unit_length; + } +#ifdef changeTrianglesOwnership + allToAll (effecPress_F, &sendCellsListReversed, &recvCellsListReversed, 1); +#endif + allToAll (effecPress_F, sendCellsList_F, recvCellsList_F, 1); +} + void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) { int numProcs, me; if (reduced_comm_id != MPI_COMM_NULL) @@ -1722,7 +1736,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep std::vector > marineBdyExtensionMap; // local map to be created by importFields importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, - stiffnessFactor_F, effecPress_F, muFriction_F, bedRougnessRC_F, + stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, temperature_F, smb_F, minThickness); import2DFieldsObservations(marineBdyExtensionMap, diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp index 7f27e9aab6cc..13cf05897f3c 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp @@ -90,7 +90,7 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F, void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, - double const* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, + double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, double* const dirichletVelocityXValue = 0, double* const dirichletVelocitYValue = 0, double* u_normal_F = 0, double* bodyForce_F = 0, double* dissipation_heat_F = 0, double* xVelocityOnCell = 0, double* yVelocityOnCell = 0, double const * deltat = 0, @@ -162,8 +162,9 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices, const std::vector& bedTopographyData, const std::vector& smbData, const std::vector& stiffnessFactorData, - const std::vector& effecPressData, + std::vector& effecPressData, const std::vector& muFrictionData, + const std::vector& bedRoughnessData, const std::vector& temperatureDataOnPrisms, std::vector& bodyForceOnBasalCell, std::vector& dissipationHeatOnPrisms, @@ -208,7 +209,7 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id); void importFields(std::vector >& marineBdyExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, - double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); + double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* bedRoughness_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); void import2DFieldsObservations(std::vector >& marineBdyExtensionMap, double const * lowerSurface_F, @@ -232,6 +233,7 @@ std::vector extendMaskByOneLayer(int const* verticesMask_F); void exportDissipationHeat(double * dissipationHeat_F); void exportBodyForce(double * bodyForce_F); void exportBeta(double * beta_F); +void exportEffectivePressure(double * effecPress_F); void get_prism_velocity_on_FEdges(double* uNormal, const std::vector& velocityOnCells, From 079213d43966c4c8703d18e15d5c4bce53cbfe45 Mon Sep 17 00:00:00 2001 From: Jeremy Brooks Date: Fri, 20 Dec 2024 12:28:30 -0800 Subject: [PATCH 3/4] Debris-friction slip law interface Adding two fields (bulkFriction, basalDebris) to use with debris-bed friction slip law. --- .../Interface_velocity_solver.cpp | 23 +++++++++++++------ .../Interface_velocity_solver.hpp | 9 +++++--- 2 files changed, 22 insertions(+), 10 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp index dcc3192f64c0..0cb56962f6b4 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp @@ -69,7 +69,7 @@ std::vector indexToTriangleID, std::vector indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge, fVertexToTriangleID, fCellToVertex, iceMarginEdgesLIds, dirichletNodesIDs; std::vector dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells, - elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell; + elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell; std::vector isVertexBoundary, isBoundaryEdge; // only needed for creating ASCII mesh @@ -241,7 +241,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, double* effecPress_F, double const* muFriction_F, - double const* bedRoughnessRC_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double* const dirichletVelocityXValue, double* const dirichletVelocitYValue, double* u_normal_F, double* bodyForce_F, double* dissipation_heat_F, double* xVelocityOnCell, double* yVelocityOnCell, double const* deltat, @@ -270,7 +270,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower if (!isDomainEmpty) { std::vector > marineBdyExtensionMap; - importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, temperature_F, smb_F, minThickness); + importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, bulkFriction_F, basalDebris_F, temperature_F, smb_F, minThickness); std::vector regulThk(thicknessData); for (int index = 0; index < nVertices; index++) @@ -287,7 +287,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta, regulThk, levelsNormalizedThickness, elevationData, thicknessData, betaData, bedTopographyData, smbData, - stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, + stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData, temperatureDataOnPrisms, bodyForceOnBasalCell, dissipationHeatOnPrisms, velocityOnVertices, albany_error, dt); *error=albany_error; @@ -1120,7 +1120,7 @@ double signedTriangleAreaOnSphere(const double* x, const double* y, void importFields(std::vector >& marineBdyExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F, double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, - double const* bedRoughnessRC_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double const * temperature_F, double const * smb_F, double eps) { int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1; @@ -1137,6 +1137,10 @@ void importFields(std::vector >& marineBdyExtensionMap, dou muFrictionData.assign(nVertices, 1e10); if (bedRoughnessRC_F!= 0) bedRoughnessData.assign(nVertices, 1e10); + if (bulkFriction_F != 0) + bulkFrictionData.assign(nVertices, 1e10); + if (basalDebris_F != 0) + basalDebrisData.assign(nVertices, 1e10); if(temperature_F != 0) temperatureDataOnPrisms.assign(nLayers * nTriangles, 1e10); if (smb_F != 0) @@ -1161,6 +1165,10 @@ void importFields(std::vector >& marineBdyExtensionMap, dou muFrictionData[index] = muFriction_F[iCell]; if (bedRoughnessRC_F != 0) bedRoughnessData[index] = bedRoughnessRC_F[iCell] / unit_length; + if (bulkFriction_F != 0) + bulkFrictionData[index] = bulkFriction_F[iCell]; + if (basalDebris_F != 0) + basalDebrisData[index] = basalDebris_F[iCell]; } int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles; @@ -1650,7 +1658,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep double const* surfaceAirTemperature_F, double const* basalHeatFlux_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, - double const* bedRoughnessRC_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double const* thickness_F, double const* thicknessUncertainty_F, double const* smb_F, double const* smbUncertainty_F, double const* bmb_F, double const* bmbUncertainty_F, @@ -1736,7 +1744,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep std::vector > marineBdyExtensionMap; // local map to be created by importFields importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, - stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, + stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, bulkFriction_F, basalDebris_F, temperature_F, smb_F, minThickness); import2DFieldsObservations(marineBdyExtensionMap, @@ -1803,6 +1811,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep write_ascii_mesh_field(effecPressData, "effective_pressure"); + // These two fields are more complicated than the others so cannot use the function to write out std::cout << "Writing temperature.ascii." << std::endl; diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp index 13cf05897f3c..31db53c1686e 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp @@ -90,7 +90,7 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F, void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, - double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, + double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double* const dirichletVelocityXValue = 0, double* const dirichletVelocitYValue = 0, double* u_normal_F = 0, double* bodyForce_F = 0, double* dissipation_heat_F = 0, double* xVelocityOnCell = 0, double* yVelocityOnCell = 0, double const * deltat = 0, @@ -134,7 +134,7 @@ void write_ascii_mesh(int const* indexToCellID_F, double const* surfaceAirTemperature_F, double const* basalHeatFlux_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, - double const* bedRoughnessRC_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double const* thickness_F, double const* thicknessUncertainty_F, double const* smb_F, double const* smbUncertainty_F, double const* bmb_F, double const* bmbUncertainty_F, @@ -165,6 +165,8 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices, std::vector& effecPressData, const std::vector& muFrictionData, const std::vector& bedRoughnessData, + const std::vector& bulkFrictionData, + const std::vector& basalDebrisData, const std::vector& temperatureDataOnPrisms, std::vector& bodyForceOnBasalCell, std::vector& dissipationHeatOnPrisms, @@ -209,7 +211,8 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id); void importFields(std::vector >& marineBdyExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, - double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* bedRoughness_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); + double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* bedRoughness_F = 0, + double const* bulkFriction_F = 0, double const* basalDebris_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); void import2DFieldsObservations(std::vector >& marineBdyExtensionMap, double const * lowerSurface_F, From 4b67208c4e0c4fa68373980af87d299ea5e22098 Mon Sep 17 00:00:00 2001 From: Jeremy Brooks Date: Mon, 18 Aug 2025 12:30:12 -0700 Subject: [PATCH 4/4] Adding spatial debris-bed friction fields to MPAS Updated Registry and MPAS velocity interface --- components/mpas-albany-landice/src/Registry.xml | 8 ++++++++ .../src/mode_forward/mpas_li_velocity_external.F | 12 ++++++++++-- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 3764768ba770..0176f9f8dfc4 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -827,6 +827,8 @@ + + @@ -1642,6 +1644,12 @@ is the value of that variable from the *previous* time level! + + diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F index b5eb319655cd..8eb852cec546 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F @@ -453,6 +453,8 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited real (kind=RKIND), dimension(:), pointer :: muFriction real (kind=RKIND), dimension(:), pointer :: bedRoughnessRC + real (kind=RKIND), dimension(:), pointer :: bulkFriction + real (kind=RKIND), dimension(:), pointer :: basalDebrisFactor real (kind=RKIND), pointer :: deltat integer, dimension(:), pointer :: vertexMask, cellMask, edgeMask integer, dimension(:,:), pointer :: dirichletVelocityMask @@ -518,6 +520,8 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(velocityPool, 'drivingStress', drivingStress) call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) call mpas_pool_get_array(velocityPool, 'bedRoughnessRC', bedRoughnessRC) + call mpas_pool_get_array(velocityPool, 'bulkFriction', bulkFriction) + call mpas_pool_get_array(velocityPool, 'basalDebrisFactor', basalDebrisFactor) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) @@ -600,7 +604,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, & betaSolve, sfcMassBal, temperature, stiffnessFactor, & effectivePressureLimited, muFriction, & - bedRoughnessRC, & + bedRoughnessRC, bulkFriction, basalDebrisFactor, & uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1 normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values deltat, albanyVelocityError) ! return values @@ -791,6 +795,8 @@ subroutine li_velocity_external_write_albany_mesh(domain) real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: muFriction real (kind=RKIND), dimension(:), pointer :: bedRoughnessRC + real (kind=RKIND), dimension(:), pointer :: bulkFriction + real (kind=RKIND), dimension(:), pointer :: basalDebrisFactor integer, dimension(:,:), pointer :: dirichletVelocityMask type (mpas_pool_type), pointer :: meshPool, geometryPool, thermalPool, observationsPool, velocityPool, scratchPool, hydroPool real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density @@ -853,6 +859,8 @@ subroutine li_velocity_external_write_albany_mesh(domain) call mpas_pool_get_array(velocityPool, 'beta', beta) call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) call mpas_pool_get_array(velocityPool, 'bedRoughnessRC', bedRoughnessRC) + call mpas_pool_get_array(velocityPool, 'bulkFriction', bulkFriction) + call mpas_pool_get_array(velocityPool, 'basalDebrisFactor', basalDebrisFactor) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) @@ -906,7 +914,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) surfaceAirTemperature, basalHeatFlux, & stiffnessFactor, & effectivePressure, muFriction, & - bedRoughnessRC, & + bedRoughnessRC, bulkFriction, basalDebrisFactor, & thickness, thicknessUncertainty, & sfcMassBal, sfcMassBalUncertainty, & floatingBasalMassBal, floatingBasalMassBalUncertainty, &