Skip to content

Add debris-bed friction slip law #135

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

Open
wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Add debris-bed friction slip law #135

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fa6b5ca
Add new spatial bed roughness variable for reg. Coulomb friction
matthewhoffman Oct 29, 2024
1c01189
Fixing typos and exporting effective pressure field
mperego Nov 1, 2024
079213d
Debris-friction slip law interface
jeremy-brooks-1 Dec 20, 2024
4b67208
Adding spatial debris-bed friction fields to MPAS
jeremy-brooks-1 Aug 18, 2025
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
12 changes: 12 additions & 0 deletions components/mpas-albany-landice/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -826,6 +826,9 @@
<!-- They do not need to be in the input file for the SIA dycore, and will be ignored if they are. -->
<var name="beta" packages="higherOrderVelocity"/>
<var name="muFriction" packages="higherOrderVelocity"/>
<var name="bedRoughnessRC" packages="higherOrderVelocity"/>
<var name="bulkFriction" packages="higherOrderVelocity"/>
<var name="basalDebrisFactor" packages="higherOrderVelocity"/>
<var name="effectivePressure" packages="higherOrderVelocity"/>
<var name="dirichletVelocityMask" packages="higherOrderVelocity"/>
<var name="uReconstructX" packages="higherOrderVelocity"/>
Expand Down Expand Up @@ -1638,6 +1641,15 @@ is the value of that variable from the *previous* time level!
<var name="stiffnessFactor" type="real" dimensions="nCells Time" default_value="1.0"
description="stiffness factor applied to effective viscosity in higher-order velocity solve"
/>
<var name="bedRoughnessRC" type="real" dimensions="nCells Time" units="m" default_value="1.0"
description="input value of bed roughness parameter for regularized Coulomb basal friction. This corresponds to lambda_max / m_max."
/>
<var name="bulkFriction" type="real" dimensions="nCells Time" default_value="0.0"
description="input value of bulk friction coefficient in debris-bed friction sliding law."
/>
<var name="basalDebrisFactor" type="real" dimensions="nCells Time" units="Pa yr m^-1" default_value="0.0"
Copy link

@matthewhoffman matthewhoffman Aug 18, 2025

Choose a reason for hiding this comment

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

@mperego , are these the correct units we should use here? Do we need to account for any scaling that Albany does here? (noting that I'm not sure if Jeremy had implemented scaling on these fields on the Albany side)

Copy link

@mperego mperego Aug 18, 2025

Choose a reason for hiding this comment

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

Based on how they are used in Albany, I think bedRoughnessRC is fine, but basalDebrisFactor should be in kPa yr m^{-1}, like muFriction. We can convert the unis in the interface_velocity_solver.cpp file.

description="input value of basal debris factor in debris-bed friction sliding law."
/>
<var name="exx" type="real" dimensions="nCells Time" units="s^-1"
description="x-component of depth-integrated strain rate"
/>
Expand Down
42 changes: 37 additions & 5 deletions components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ std::vector<int> indexToTriangleID,
std::vector<int> indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge,
fVertexToTriangleID, fCellToVertex, iceMarginEdgesLIds, dirichletNodesIDs;
std::vector<double> dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells,
elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell;
elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell;
std::vector<bool> isVertexBoundary, isBoundaryEdge;

// only needed for creating ASCII mesh
Expand Down Expand Up @@ -240,7 +240,8 @@ 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* 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,
Expand Down Expand Up @@ -269,7 +270,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower
if (!isDomainEmpty) {

std::vector<std::pair<int, int> > 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, bedRoughnessRC_F, bulkFriction_F, basalDebris_F, temperature_F, smb_F, minThickness);

std::vector<double> regulThk(thicknessData);
for (int index = 0; index < nVertices; index++)
Expand All @@ -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, bulkFrictionData, basalDebrisData,
temperatureDataOnPrisms, bodyForceOnBasalCell, dissipationHeatOnPrisms, velocityOnVertices,
albany_error, dt);
*error=albany_error;
Expand All @@ -299,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);

Expand Down Expand Up @@ -1117,6 +1120,7 @@ double signedTriangleAreaOnSphere(const double* x, const double* y,

void importFields(std::vector<std::pair<int, int> >& 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* bulkFriction_F, double const* basalDebris_F,
double const * temperature_F, double const * smb_F, double eps) {

int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1;
Expand All @@ -1131,6 +1135,12 @@ void importFields(std::vector<std::pair<int, int> >& marineBdyExtensionMap, dou
effecPressData.assign(nVertices, 1e10);
if (muFriction_F!= 0)
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)
Expand All @@ -1153,6 +1163,12 @@ void importFields(std::vector<std::pair<int, int> >& 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] / 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;
Expand Down Expand Up @@ -1361,6 +1377,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)
Expand Down Expand Up @@ -1630,6 +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* 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,
Expand Down Expand Up @@ -1715,7 +1744,8 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep

std::vector<std::pair<int, int> > 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, bedRoughnessRC_F, bulkFriction_F, basalDebris_F,
temperature_F, smb_F, minThickness);

import2DFieldsObservations(marineBdyExtensionMap,
thicknessUncertainty_F,
Expand Down Expand Up @@ -1760,6 +1790,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");
Expand All @@ -1780,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;
Expand Down
12 changes: 9 additions & 3 deletions components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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* 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,
Expand Down Expand Up @@ -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* 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,
Expand Down Expand Up @@ -161,8 +162,11 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices,
const std::vector<double>& bedTopographyData,
const std::vector<double>& smbData,
const std::vector<double>& stiffnessFactorData,
const std::vector<double>& effecPressData,
std::vector<double>& effecPressData,
Copy link

@matthewhoffman matthewhoffman Aug 18, 2025

Choose a reason for hiding this comment

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

@mperego , was this change intentional?

Copy link

@mperego mperego Aug 18, 2025

Choose a reason for hiding this comment

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

Yes. Now Albany can return the effective pressure used (e.g., the open-connection one).

const std::vector<double>& muFrictionData,
const std::vector<double>& bedRoughnessData,
const std::vector<double>& bulkFrictionData,
const std::vector<double>& basalDebrisData,
const std::vector<double>& temperatureDataOnPrisms,
std::vector<double>& bodyForceOnBasalCell,
std::vector<double>& dissipationHeatOnPrisms,
Expand Down Expand Up @@ -207,7 +211,8 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id);

void importFields(std::vector<std::pair<int, int> >& 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* 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<std::pair<int, int> >& marineBdyExtensionMap,
double const * lowerSurface_F,
Expand All @@ -231,6 +236,7 @@ std::vector<int> 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<double>& velocityOnCells,
Expand Down
14 changes: 14 additions & 0 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -452,6 +452,9 @@ 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), dimension(:), pointer :: bulkFriction
real (kind=RKIND), dimension(:), pointer :: basalDebrisFactor
real (kind=RKIND), pointer :: deltat
integer, dimension(:), pointer :: vertexMask, cellMask, edgeMask
integer, dimension(:,:), pointer :: dirichletVelocityMask
Expand Down Expand Up @@ -516,6 +519,9 @@ 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, '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)
Expand Down Expand Up @@ -598,6 +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, bulkFriction, basalDebrisFactor, &
uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1
normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values
deltat, albanyVelocityError) ! return values
Expand Down Expand Up @@ -787,6 +794,9 @@ 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
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
Expand Down Expand Up @@ -848,6 +858,9 @@ 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, '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)

Expand Down Expand Up @@ -901,6 +914,7 @@ subroutine li_velocity_external_write_albany_mesh(domain)
surfaceAirTemperature, basalHeatFlux, &
stiffnessFactor, &
effectivePressure, muFriction, &
bedRoughnessRC, bulkFriction, basalDebrisFactor, &
thickness, thicknessUncertainty, &
sfcMassBal, sfcMassBalUncertainty, &
floatingBasalMassBal, floatingBasalMassBalUncertainty, &
Expand Down