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 1 commit
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
1 change: 1 addition & 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,7 @@
<!-- 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="effectivePressure" packages="higherOrderVelocity"/>
<var name="dirichletVelocityMask" packages="higherOrderVelocity"/>
<var name="uReconstructX" packages="higherOrderVelocity"/>
Expand Down
24 changes: 19 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 @@ -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,
Expand Down Expand Up @@ -270,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, 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<double> regulThk(thicknessData);
for (int index = 0; index < nVertices; index++)
Expand All @@ -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;
Expand All @@ -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);

Expand Down Expand Up @@ -1158,7 +1160,7 @@ void importFields(std::vector<std::pair<int, int> >& 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;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -1722,7 +1736,7 @@ 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, bedRougnessRC_F,
stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F,
temperature_F, smb_F, minThickness);

import2DFieldsObservations(marineBdyExtensionMap,
Expand Down
8 changes: 5 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 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,
Expand Down Expand Up @@ -162,8 +162,9 @@ 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>& temperatureDataOnPrisms,
std::vector<double>& bodyForceOnBasalCell,
std::vector<double>& dissipationHeatOnPrisms,
Expand Down Expand Up @@ -208,7 +209,7 @@ 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* 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 @@ -232,6 +233,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