Skip to content

Changed platform erosion to take cell arg and started parallelisation #171

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
wilfchun wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Changed platform erosion to take cell arg and started parallelisation #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
61 changes: 39 additions & 22 deletions src/do_shore_platform_erosion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ using std::array;
using std::shuffle;

#include "cme.h"
#include "cell.h"
#include "hermite_cubic.h"
#include "simulation.h"
#include "coast.h"
Expand Down Expand Up @@ -138,13 +139,20 @@ int CSimulation::nDoAllShorePlatFormErosion(void)
FillInBeachProtectionHolesAndRemoveLegacyCliffs();

// Finally calculate actual platform erosion on all sea cells (both on profiles, and between profiles)
#pragma omp parallel for schedule(dynamic) collapse(2) \
reduction(+:m_dThisIterActualPlatformErosionFineCons, m_dThisIterFineSedimentToSuspension, \
m_dThisIterActualPlatformErosionSandCons, m_dThisIterActualPlatformErosionCoarseCons, \
m_ulThisIterNumActualPlatformErosionCells, m_dDepositionSandDiff, m_dDepositionCoarseDiff)
for (int nX = 0; nX < m_nXGridSize; nX++)
{
for (int nY = 0; nY < m_nYGridSize; nY++)
{
if (m_pRasterGrid->m_Cell[nX][nY].bPotentialPlatformErosion())
// Cache cell reference to avoid repeated array lookups
auto& rCell = m_pRasterGrid->m_Cell[nX][nY];

if (rCell.bPotentialPlatformErosion())
// Calculate actual (supply-limited) shore platform erosion on each cell that has potential platform erosion, also add the eroded sand/coarse sediment to that cell's polygon, ready to be redistributed within the polygon during beach erosion/deposition
DoActualPlatformErosionOnCell(nX, nY);
DoActualPlatformErosionOnCell(nX, nY, rCell);
}
}

Expand Down Expand Up @@ -732,22 +740,22 @@ int CSimulation::nCalcPotentialPlatformErosionBetweenProfiles(int const nCoast,
//===============================================================================================================================
//! Calculates actual (constrained by available sediment) erosion of the consolidated shore platform on a single sea cell
//===============================================================================================================================
void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY, CGeomCell& rCell)
{
// LogStream << m_ulIter << ": doing platform erosion on cell [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;

// Get the beach protection factor, which quantifies the extent to which unconsolidated sediment on the shore platform (beach) protects the shore platform
double const dBeachProtectionFactor = m_pRasterGrid->m_Cell[nX][nY].dGetBeachProtectionFactor();
double const dBeachProtectionFactor = rCell.dGetBeachProtectionFactor();

if (bFPIsEqual(dBeachProtectionFactor, 0.0, TOLERANCE))
// The beach is sufficiently thick to prevent any platform erosion (or we are down to basement)
return;

// Get the potential depth of potential erosion, considering beach protection
double const dThisPotentialErosion = m_pRasterGrid->m_Cell[nX][nY].dGetPotentialPlatformErosion() * dBeachProtectionFactor;
double const dThisPotentialErosion = rCell.dGetPotentialPlatformErosion() * dBeachProtectionFactor;

// We will be eroding the topmost layer that has non-zero thickness
int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
int const nThisLayer = rCell.nGetTopNonZeroLayerAboveBasement();

// Safety check
if (nThisLayer == INT_NODATA)
Expand All @@ -761,9 +769,9 @@ void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
return;

// OK, we have a layer that can be eroded so find out how much consolidated sediment we have available on this cell
double const dExistingAvailableFine = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetFineDepth();
double const dExistingAvailableSand = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetSandDepth();
double const dExistingAvailableCoarse = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetCoarseDepth();
double const dExistingAvailableFine = rCell.pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetFineDepth();
double const dExistingAvailableSand = rCell.pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetSandDepth();
double const dExistingAvailableCoarse = rCell.pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetCoarseDepth();

// Now partition the total lowering for this cell between the three size fractions: do this by relative erodibility
int const nFineWeight = (dExistingAvailableFine > 0 ? 1 : 0);
Expand All @@ -787,7 +795,7 @@ void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
dTotActualErosion += dFineEroded;

// Set the value for this layer
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetFineDepth(dRemaining);
rCell.pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetFineDepth(dRemaining);

// And set the changed-this-timestep switch
m_bConsChangedThisIter[nThisLayer] = true;
Expand All @@ -809,10 +817,10 @@ void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
dTotActualErosion += dSandEroded;

// Set the new value of sand consolidated sediment depth for this layer
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetSandDepth(dRemaining);
rCell.pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetSandDepth(dRemaining);

// And add this to the depth of sand unconsolidated sediment for this layer
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->AddSandDepth(dSandEroded);
rCell.pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->AddSandDepth(dSandEroded);

// Set the changed-this-timestep switch
m_bConsChangedThisIter[nThisLayer] = true;
Expand All @@ -833,10 +841,10 @@ void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
dTotActualErosion += dCoarseEroded;

// Set the new value of coarse consolidated sediment depth for this layer
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetCoarseDepth(dRemaining);
rCell.pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetCoarseDepth(dRemaining);

// And add this to the depth of coarse unconsolidated sediment for this layer
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->AddCoarseDepth(dCoarseEroded);
rCell.pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->AddCoarseDepth(dCoarseEroded);

// Set the changed-this-timestep switch
m_bConsChangedThisIter[nThisLayer] = true;
Expand All @@ -849,26 +857,31 @@ void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
if (dTotActualErosion > 0)
{
// We did, so set the actual erosion value for this cell
m_pRasterGrid->m_Cell[nX][nY].SetActualPlatformErosion(dTotActualErosion);
rCell.SetActualPlatformErosion(dTotActualErosion);

// Recalculate the elevation of every layer
m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
rCell.CalcAllLayerElevsAndD50();

// And update the cell's sea depth
m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
rCell.SetSeaDepth();

// Update per-timestep totals
m_ulThisIterNumActualPlatformErosionCells++;

// Add eroded sand/coarse sediment for this cell to the polygon that contains the cell, ready for redistribution during beach erosion/deposition (fine sediment has already been dealt with)
int nPolyID = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
int nPolyCoastID = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonCoastID();
int nPolyID = rCell.nGetPolygonID();
int nPolyCoastID = rCell.nGetPolygonCoastID();

if (nPolyID == INT_NODATA)
{
// Can get occasional problems with polygon rasterization near the coastline, so also search the eight adjacent cells
array<int, 8> nDirection = {NORTH, NORTH_EAST, EAST, SOUTH_EAST, SOUTH, SOUTH_WEST, WEST, NORTH_WEST};
shuffle(nDirection.begin(), nDirection.end(), m_Rand[0]);

// Thread-safe: protect RNG access in parallel region
#pragma omp critical(rng_shuffle)
{
shuffle(nDirection.begin(), nDirection.end(), m_Rand[0]);
}

for (int n = 0; n < 8; n++)
{
Expand Down Expand Up @@ -1023,8 +1036,12 @@ void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
}

// All OK, so add this to the polygon's total of unconsolidated sand/coarse sediment, to be deposited or moved later. These values are +ve (deposition)
m_VCoast[nPolyCoastID].pGetPolygon(nPolyID)->AddPlatformErosionUnconsSand(dSandEroded);
m_VCoast[nPolyCoastID].pGetPolygon(nPolyID)->AddPlatformErosionUnconsCoarse(dCoarseEroded);
// Thread-safe: protect polygon updates in parallel region
#pragma omp critical(polygon_updates)
{
m_VCoast[nPolyCoastID].pGetPolygon(nPolyID)->AddPlatformErosionUnconsSand(dSandEroded);
m_VCoast[nPolyCoastID].pGetPolygon(nPolyID)->AddPlatformErosionUnconsCoarse(dCoarseEroded);
}
}
}

Expand Down
3 changes: 2 additions & 1 deletion src/simulation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ class CGeomCoastPolygon;
class CRWCliff;
class CRWSedInputEvent;
class CRWCellLandform;
class CGeomCell;

class CSimulation
{
Expand Down Expand Up @@ -1694,7 +1695,7 @@ class CSimulation
double dCalcBeachProtectionFactor(int const, int const, double const);
void FillInBeachProtectionHolesAndRemoveLegacyCliffs(void);
void FillPotentialPlatformErosionHoles(void);
void DoActualPlatformErosionOnCell(int const, int const);
void DoActualPlatformErosionOnCell(int const, int const, CGeomCell&);
double dLookUpErosionPotential(double const);
static CGeom2DPoint PtChooseEndPoint(int const, CGeom2DPoint const*, CGeom2DPoint const*, double const, double const, double const, double const);
int nGetCoastNormalEndPoint(int const, int const, int const, CGeom2DPoint const*, double const, CGeom2DPoint *, CGeom2DIPoint *, bool const);
Expand Down