CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
do_shore_platform_erosion.cpp
Go to the documentation of this file.
1
12
13/* ==============================================================================================================================
14
15 This file is part of CoastalME, the Coastal Modelling Environment.
16
17 CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
20
21 You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22
23==============================================================================================================================*/
24#include <assert.h>
25
26#include <cmath>
27
28#include <iostream>
29using std::cerr;
30using std::endl;
31using std::ios;
32
33#include <array>
34using std::array;
35
36#include <algorithm>
37using std::shuffle;
38
39#include "cme.h"
40#include "hermite_cubic.h"
41#include "simulation.h"
42#include "coast.h"
43#include "2di_point.h"
44
45//===============================================================================================================================
47//===============================================================================================================================
49{
51 LogStream << m_ulIter << ": Calculating shore platform erosion" << endl;
52
53 // // DEBUG CODE ===========================================================================================================
54 // string strOutFile = m_strOutPath + "01_polygon_shore_platform_test_";
55 // strOutFile += to_string(m_ulIter);
56 // strOutFile += ".tif";
57 //
58 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
59 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
60 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
61 // pDataSet->SetGeoTransform(m_dGeoTransform);
62 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
63 // int n = 0;
64 // int nInPoly = 0;
65 // int nNotInPoly = 0;
66 //
67 // for (int nY = 0; nY < m_nYGridSize; nY++)
68 // {
69 // for (int nX = 0; nX < m_nXGridSize; nX++)
70 // {
71 // int nID = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
72 // if (nID == INT_NODATA)
73 // nNotInPoly++;
74 // else
75 // nInPoly++;
76 //
77 // pdRaster[n++] = nID;
78 // }
79 // }
80 //
81 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
82 // pBand->SetNoDataValue(m_dMissingValue);
83 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
84 // if (nRet == CE_Failure)
85 // LogStream << nRet << endl;
86 //
87 // GDALClose(pDataSet);
88 // delete[] pdRaster;
89 //
90 // LogStream << m_ulIter << " Number of cells in a polygon = " << nInPoly << endl;
91 // LogStream << m_ulIter << " Number of cells not in any polygon = " << nNotInPoly << endl;
92 // // DEBUG CODE ===========================================================================================================
93
94 // TODO 023 Only do potential erosion if cell is in a polygon
95
96 // Set direction
97 static bool bDownCoast = true;
98
99 // Do this for each coast
100 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
101 {
102 int const nNumProfiles = m_VCoast[nCoast].nGetNumProfiles();
103
104 // Calculate potential platform erosion along each coastline-normal profile, do in down-coast sequence
105 for (int nn = 0; nn < nNumProfiles; nn++)
106 {
107 CGeomProfile *pProfile;
108
109 if (bDownCoast)
110 pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(nn);
111
112 else
113 pProfile = m_VCoast[nCoast].pGetProfileWithUpCoastSeq(nn);
114
115 int const nRet = nCalcPotentialPlatformErosionOnProfile(nCoast, pProfile);
116
117 if (nRet != RTN_OK)
118 return nRet;
119 }
120
121 // Calculate potential platform erosion between the coastline-normal profiles. Do this in down-coast sequence
122 for (int nn = 0; nn < nNumProfiles - 1; nn++)
123 {
124 CGeomProfile *pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(nn);
125
126 // Calculate potential erosion for sea cells between this profile and the next profile (or up to the edge of the grid) on these cells
128
129 if (nRet != RTN_OK)
130 return nRet;
131
133
134 if (nRet != RTN_OK)
135 return nRet;
136 }
137 }
138
139 // Swap direction for next timestep
140 bDownCoast = ! bDownCoast;
141
142 // Fills in 'holes' in the potential platform erosion i.e. orphan cells which get omitted because of rounding problems
144
145 // Do the same for beach protection
147
148 // Finally calculate actual platform erosion on all sea cells (both on profiles, and between profiles)
149 for (int nX = 0; nX < m_nXGridSize; nX++)
150 {
151 for (int nY = 0; nY < m_nYGridSize; nY++)
152 {
153 if (m_pRasterGrid->m_Cell[nX][nY].bPotentialPlatformErosion())
154 // 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
156 }
157 }
158
160 {
161 LogStream << m_ulIter << ": total potential shore platform erosion (m^3) = " << m_dThisIterPotentialPlatformErosion * m_dCellArea << " (on profiles = " << m_dTotPotentialPlatformErosionOnProfiles * m_dCellArea << ", between profiles = " << m_dTotPotentialPlatformErosionBetweenProfiles * m_dCellArea << ")" << endl;
162
164 }
165
166 return RTN_OK;
167}
168
169//===============================================================================================================================
171//===============================================================================================================================
173{
174 // Only work on this profile if it is problem-free TODO 024 Or if it has just hit dry land?
175 if (! pProfile->bOKIncStartAndEndOfCoast()) // || (pProfile->nGetProblemCode() == PROFILE_DRYLAND))
176 return RTN_OK;
177
178 // Get the length of the profile (in cells) and the index of the coast point at which this profile starts
179 int const nProfSize = pProfile->nGetNumCellsInProfile();
180 int const nCoastPoint = pProfile->nGetCoastPoint();
181
182 // Get the breaking depth for this profile from the coastline point
183 double const dDepthOfBreaking = m_VCoast[nCoast].dGetDepthOfBreaking(nCoastPoint);
184
185 if (bFPIsEqual(dDepthOfBreaking, DBL_NODATA, TOLERANCE))
186 // This profile is not in the active zone, so no platform erosion here
187 return RTN_OK;
188
189 if (bFPIsEqual(dDepthOfBreaking, 0.0, TOLERANCE))
190 {
191 // Safety check, altho' this shouldn't happen
193 LogStream << m_ulIter << ": depth of breaking is zero for profile " << pProfile->nGetProfileID() << " of coast " << nCoast << endl;
194
195 return RTN_OK;
196 }
197
198 // LogStream << m_ulIter << ": ON PROFILE nProfile = " << nProfile << " dDepthOfBreaking = " << dDepthOfBreaking << endl;
199
200 // Get the height of the associated breaking wave from the coast point: this height is used in beach protection calcs
201 double const dBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nCoastPoint);
202
203 // Calculate the length of the profile in external CRS units
204 int const nSegments = pProfile->nGetProfileSize() - 1;
205 double dProfileLenXY = 0;
206
207 for (int nSeg = 0; nSeg < nSegments; nSeg++)
208 {
209 // Do once for every line segment
210 double const dSegStartX = pProfile->pPtGetPointInProfile(nSeg)->dGetX();
211 double const dSegStartY = pProfile->pPtGetPointInProfile(nSeg)->dGetY();
212 double const dSegEndX = pProfile->pPtGetPointInProfile(nSeg + 1)->dGetX(); // Is OK
213 double const dSegEndY = pProfile->pPtGetPointInProfile(nSeg + 1)->dGetY();
214
215 double const dSegLen = hypot(dSegStartX - dSegEndX, dSegStartY - dSegEndY);
216 dProfileLenXY += dSegLen;
217 }
218
219 // Next calculate the average distance between profile points, again in external CRS units. Assume that the sample points are equally spaced along the profile (not quite true)
220 double const dSpacingXY = dProfileLenXY / (nProfSize - 1);
221
222 // Set up vectors for the coastline-normal profile elevations. The length of this vector line is given by the number of cells 'under' the profile. Thus each point on the vector relates to a single cell in the grid. This assumes that all points on the profile vector are equally spaced (not quite true, depends on the orientation of the line segments which comprise the profile)
223 // The elevation of each of these profile points is the elevation of the centroid of the cell that is 'under' the point. However we cannot always be confident that this is the 'true' elevation of the point on the vector since (unless the profile runs planview N-S or W-E) the vector does not always run exactly through the centroid of the cell
224 vector<double> VdProfileZ(nProfSize, 0); // Initial (pre-erosion) elevation of both consolidated and unconsolidated sediment for cells 'under' the profile
225 vector<double> VdProfileDistXY(nProfSize, 0); // Along-profile distance measured from the coast, in external CRS units
226 vector<double> dVConsProfileZ(nProfSize, 0); // Initial (pre-erosion) elevation of consolidated sediment only for cells 'under' the profile
227 vector<double> dVConsZDiff(nProfSize, 0);
228 vector<double> dVConsSlope(nProfSize, 0);
229
230 for (int i = 0; i < nProfSize; i++)
231 {
232 int const nX = pProfile->pPtiVGetCellsInProfile()->at(i).nGetX();
233 int const nY = pProfile->pPtiVGetCellsInProfile()->at(i).nGetY();
234
235 // Get the number of the highest layer with non-zero thickness
236 int const nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
237
238 // Safety check
239 if (nTopLayer == INT_NODATA)
241
242 if (nTopLayer == NO_NONZERO_THICKNESS_LAYERS)
243 // TODO 025 We are down to basement
244 return RTN_OK;
245
246 // Get the elevation for consolidated sediment only on this cell
247 dVConsProfileZ[i] = m_pRasterGrid->m_Cell[nX][nY].dGetConsSedTopForLayerAboveBasement(nTopLayer);
248
249 // Get the elevation for both consolidated and unconsolidated sediment on this cell
250 VdProfileZ[i] = m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
251
252 // And store the X-Y plane distance from the start of the profile
253 VdProfileDistXY[i] = i * dSpacingXY;
254 }
255
256 for (int i = 0; i < nProfSize - 1; i++)
257 {
258 // For the consolidated-only profile, get the Z differences (already in external CRS units)
259 dVConsZDiff[i] = dVConsProfileZ[i] - dVConsProfileZ[i + 1];
260
261 // Calculate dZ/dXY, the Z slope (i.e. rise over run) in the XY direction. Note that we use the elevation difference on the seaward side of 'this' point
262 dVConsSlope[i] = dVConsZDiff[i] / dSpacingXY;
263 }
264
265 // Sort out the final value
266 dVConsSlope[nProfSize - 1] = dVConsSlope[nProfSize - 2];
267
269 {
270 // Smooth the vector of slopes for the consolidated-only profile
271 dVConsSlope = dVSmoothProfileSlope(&dVConsSlope);
272 }
273
274 vector<double> dVProfileDepthOverDB(nProfSize, 0); // Depth over wave breaking depth at the coastline-normal sample points
275 vector<double> dVProfileErosionPotential(nProfSize, 0); // Erosion potential at the coastline-normal sample points
276
277 // Calculate the erosion potential along this profile using the shape function
278 double dTotalErosionPotential = 0;
279
280 for (int i = 0; i < nProfSize; i++)
281 {
282 // Use the actual depth of water here (i.e. the depth to the top of the unconsolidated sediment, including the thickness of consolidated sediment beneath it)
283 dVProfileDepthOverDB[i] = m_dThisIterSWL - VdProfileZ[i];
284 dVProfileDepthOverDB[i] /= dDepthOfBreaking;
285
286 // Constrain dDepthOverDB[i] to be between 0 (can get small -ve values due to rounding errors) and m_dDepthOverDBMax
287 dVProfileDepthOverDB[i] = tMax(dVProfileDepthOverDB[i], 0.0);
288 dVProfileDepthOverDB[i] = tMin(dVProfileDepthOverDB[i], m_dDepthOverDBMax);
289
290 // And then use the look-up table to find the value of erosion potential at this point on the profile
291 dVProfileErosionPotential[i] = dLookUpErosionPotential(dVProfileDepthOverDB[i]);
292
293 // If erosion potential (a -ve value) is tiny, set it to zero
294 if (dVProfileErosionPotential[i] > -SEDIMENT_ELEV_TOLERANCE)
295 dVProfileErosionPotential[i] = 0;
296
297 // Keep track of the total erosion potential for this profile
298 dTotalErosionPotential += dVProfileErosionPotential[i];
299 }
300
301 // Constrain erosion potential at every point on the profile, so that the integral of erosion potential on the whole profile is unity (Walkden and Hall 2005). Note that here, erosion potential is -ve so we must constrain to -1
302 for (int i = 0; i < nProfSize; i++)
303 {
304 if (dTotalErosionPotential < 0)
305 dVProfileErosionPotential[i] /= (-dTotalErosionPotential);
306 }
307
308 vector<double> dVRecessionXY(nProfSize, 0);
309 vector<double> dVSCAPEXY(nProfSize, 0);
310
311 // Calculate recession at every point on the coastline-normal profile
312 for (int i = 0; i < nProfSize; i++)
313 {
314 // dRecession = dForce * (dBeachProtection / dR) * dErosionPotential * dSlope * dTime where:
315 // dVRecession [m] is the landward migration distance defined in the profile relative (XY) CRS
316 // dForce is given by Equation 4 in Walkden & Hall, 2005
317 // dVBeachProtection [1] is beach protection factor [1, 0] = [no protection, fully protected]. (This is calculated later, see dCalcBeachProtectionFactor())
318 // dVR [m^(9/4)s^(2/3)] is the material strength and some hydrodynamic constant
319 // dVProfileErosionPotential [?] is the erosion potential at each point along the profile
320 // dVSlope [1] is the along-profile slope
321 // m_dTimeStep * 3600 [s] is the time interval in seconds
322 //
323 // dRecession is horizontal recession (along the XY direction):
324 //
325 // dVRecessionXY[i] = (dForce * dVBeachProtection[i] * dVErosionPotentialFunc[i] * dVSlope[i] * m_dTimeStep * 3600) / dVR[i]
326 //
327 // XY recession must be -ve or zero. If it is +ve then it represents accretion not erosion, which must be described by a different set of equations. So we also need to constrain XY recession to be <= 0
328 dVRecessionXY[i] = tMin(m_VCoast[nCoast].dGetWaveEnergyAtBreaking(nCoastPoint) * dVProfileErosionPotential[i] * dVConsSlope[i] / m_dR, 0.0);
329 dVSCAPEXY[i] = VdProfileDistXY[i] - dVRecessionXY[i];
330 }
331
332 vector<double> dVChangeElevZ(nProfSize, 0);
333
334 // We have calculated the XY-plane recession at every point on the profile, so now convert this to a change in Z-plane elevation at every inundated point on the profile (not the coast point). Again we use the elevation difference on the seaward side of 'this' point
335 for (int i = 1; i < nProfSize - 1; i++)
336 {
337 // Vertical lowering dZ = dXY * tan(a), where tan(a) is the slope of the SCAPE profile in the XY direction
338 double const dSCAPEHorizDist = dVSCAPEXY[i + 1] - dVSCAPEXY[i];
339 double const dSCAPEVertDist = dVConsProfileZ[i] - dVConsProfileZ[i + 1];
340 double const dSCAPESlope = dSCAPEVertDist / dSCAPEHorizDist;
341 double dDeltaZ = dVRecessionXY[i] * dSCAPESlope;
342
343 // Safety check: if thickness model has some jumps, dVConsProfileZ might be very high, limiting dSCAPESlope to 0 because all time erode a high fix quantity
344 if (dSCAPESlope > 1)
345 {
346 dDeltaZ = 0;
347 }
348
349 int const nX = pProfile->pPtiVGetCellsInProfile()->at(i).nGetX();
350 int const nY = pProfile->pPtiVGetCellsInProfile()->at(i).nGetY();
351
352 // Store the local slope of the consolidated sediment, this is just for output display purposes
353 m_pRasterGrid->m_Cell[nX][nY].SetLocalConsSlope(dVConsSlope[i]);
354
355 // dDeltaZ is zero or -ve: if dDeltaZ is zero then do nothing, if -ve then remove some sediment from this cell
356 if (dDeltaZ < 0)
357 {
358 // If there has already been potential erosion on this cell, then it must be a shared line segment (i.e. has co-incident profiles)
359 double const dPrevPotentialErosion = -m_pRasterGrid->m_Cell[nX][nY].dGetPotentialPlatformErosion();
360
361 if (dPrevPotentialErosion < 0)
362 {
363 // Average the two values
364 // LogStream << m_ulIter << ": [" << nX << "][" << nY << "] under profile " << nProfile << " has previous potential platform erosion = " << dPrevPotentialErosion << endl;
365 dDeltaZ = ((dDeltaZ + dPrevPotentialErosion) / 2);
366 }
367
368 // Constrain the lowering so we don't get negative slopes or +ve erosion amounts (dDeltaZ must be -ve), this is implicit in SCAPE
369 dDeltaZ = tMax(dDeltaZ, -dVConsZDiff[i]);
370 dDeltaZ = tMin(dDeltaZ, 0.0);
371 dVChangeElevZ[i] = dDeltaZ;
372
373 // Set the potential (unconstrained) erosion for this cell, is a +ve value
374 m_pRasterGrid->m_Cell[nX][nY].SetPotentialPlatformErosion(-dDeltaZ);
375
376 // Update this-timestep totals
378 m_dThisIterPotentialPlatformErosion -= dDeltaZ; // Since dDeltaZ is a -ve value
379 // assert(isfinite(m_dThisIterPotentialPlatformErosion));
380 // assert(m_dThisIterPotentialPlatformErosion >= 0);
381
382 // Increment the check values
385 }
386
387 // Finally, calculate the beach protection factor, this will be used in estimating actual (supply-limited) erosion
388 double const dBeachProtectionFactor = dCalcBeachProtectionFactor(nX, nY, dBreakingWaveHeight);
389 m_pRasterGrid->m_Cell[nX][nY].SetBeachProtectionFactor(dBeachProtectionFactor);
390 }
391
392 // If desired, save this coastline-normal profile data for checking purposes
394 {
395 int const nRet = nSaveProfile(nCoast, pProfile, nProfSize, &VdProfileDistXY, &dVConsProfileZ, &dVProfileDepthOverDB, &dVProfileErosionPotential, &dVConsSlope, &dVRecessionXY, &dVChangeElevZ, pProfile->pPtiVGetCellsInProfile(), &dVSCAPEXY);
396
397 if (nRet != RTN_OK)
398 return nRet;
399 }
400
401 return RTN_OK;
402}
403
404//===============================================================================================================================
406//===============================================================================================================================
407int CSimulation::nCalcPotentialPlatformErosionBetweenProfiles(int const nCoast, CGeomProfile *pProfile, int const nDirection)
408{
409 // Only work on this profile if it is problem-free
410 if (! pProfile->bOKIncStartAndEndOfCoast())
411 return RTN_OK;
412
413 int const nProfSize = pProfile->nGetNumCellsInProfile();
414 int const nCoastProfileStart = pProfile->nGetCoastPoint();
415 int const nProfileStartX = pProfile->pPtiVGetCellsInProfile()->at(0).nGetX();
416 int const nProfileStartY = pProfile->pPtiVGetCellsInProfile()->at(0).nGetY();
417 int const nCoastMax = m_VCoast[nCoast].nGetCoastlineSize();
418 int nDistFromProfile = 0;
419 int nParCoastXLast = nProfileStartX;
420 int nParCoastYLast = nProfileStartY;
421
422 // Start at the coast end of this coastline-normal profile, then move one cell forward along the coast, then construct a parallel profile from this new coastline start cell. Calculate erosion along this parallel profile in the same way as above. Move another cell forward along the coastline, do the same. Keep going until hit another profile
423 while (true)
424 {
425 // Increment the distance from the start coast-normal profile
426 nDistFromProfile++;
427
428 // Start the parallel profile nDistFromProfile cells along the coastline from the coastline-normal profile, direction depending on nDirection
429 int nThisPointOnCoast = nCoastProfileStart;
430
431 if (nDirection == DIRECTION_DOWNCOAST)
432 nThisPointOnCoast += nDistFromProfile;
433
434 else
435 nThisPointOnCoast -= nDistFromProfile;
436
437 // Have we reached the beginning of the coast?
438 if ((nDirection == DIRECTION_UPCOAST) && (nThisPointOnCoast < 0))
439 {
440 // LogStream << m_ulIter << ": LEAVING LOOP since hit nThisPointOnCoast = " << nThisPointOnCoast << " while doing potential platform erosion " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast from profile = " << nProfile << ", dist from profile = " << nDistFromProfile << endl;
441
442 break;
443 }
444
445 // Have we reached the end of the coast?
446 if ((nDirection == DIRECTION_DOWNCOAST) && (nThisPointOnCoast >= nCoastMax))
447 {
448 // LogStream << m_ulIter << ": LEAVING LOOP since hit nThisPointOnCoast = " << nThisPointOnCoast << " while doing potential platform erosion " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast from profile = " << nProfile << ", dist from profile = " << nDistFromProfile << endl;
449
450 break;
451 }
452
453 // LogStream << m_ulIter << ": from profile " << nProfile << ", doing potential platform erosion " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast, dist from profile = " << nDistFromProfile << endl;
454
455 double dDepthOfBreaking = m_VCoast[nCoast].dGetDepthOfBreaking(nThisPointOnCoast);
456
457 if (bFPIsEqual(dDepthOfBreaking, DBL_NODATA, TOLERANCE))
458 {
459 // This parallel profile is not in the active zone, so no platform erosion here. Move on to the next point along the coastline in this direction
460 // if (m_nLogFileDetail == LOG_FILE_ALL)
461 // LogStream << m_ulIter << ": not in active zone at coastline " << nCoast << " coast point " << nThisPointOnCoast << " when constructing parallel profile for potential platform erosion. Working from profile " << pProfile->nGetProfileID() << ", " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast, dist from profile = " << nDistFromProfile << endl;
462
463 continue;
464 }
465
466 // LogStream << m_ulIter << ": BETWEEN PROFILES " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast from profile " << nProfile << " nThisPointOnCoast = " << nThisPointOnCoast << " dDepthOfBreaking = " << dDepthOfBreaking << " nParProfSize = " << nParProfSize << endl;
467
468 // All is OK, so get the grid coordinates of this point, which is the coastline start point for the parallel profile
469 int const nParCoastX = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nThisPointOnCoast)->nGetX();
470 int const nParCoastY = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nThisPointOnCoast)->nGetY();
471
472 if ((nParCoastX == nParCoastXLast) && (nParCoastY == nParCoastYLast))
473 {
474 // Should not happen, but could do due to rounding errors
475 // if (m_nLogFileDetail >= LOG_FILE_ALL)
476 // LogStream << WARN << m_ulIter << ": rounding problem on coast " << nCoast << " profile " << pProfile->nGetProfileID() << " at [" << nParCoastX << "][" << nParCoastY << "]" << endl;
477
478 // So move on to the next point along the coastline in this direction
479 continue;
480 }
481
482 // Is this coastline start point the start point of an adjacent coastline-normal vector?
483 if (m_pRasterGrid->m_Cell[nParCoastX][nParCoastY].bIsProfile())
484 {
485 // LogStream << m_ulIter << ": coast " << nCoast << ", LEAVING LOOP since hit another profile at nThisPointOnCoast = " << nThisPointOnCoast << " while doing potential platform erosion " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast from profile = " << nProfile << ", dist from profile = " << nDistFromProfile << endl;
486 break;
487 }
488
489 // Get the height of the associated breaking wave from the coast point: this height is used in beach protection calcs. Note that it will be DBL_NODATA if not in active zone
490 double const dBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nThisPointOnCoast);
491
492 // OK, now construct a parallel profile
493 vector<CGeom2DIPoint> PtiVGridParProfile; // Integer coords (grid CRS) of cells under the parallel profile
494 vector<CGeom2DPoint> PtVExtCRSParProfile; // coordinates (external CRS) of cells under the parallel profile
495
496 ConstructParallelProfile(nProfileStartX, nProfileStartY, nParCoastX, nParCoastY, nProfSize, pProfile->pPtiVGetCellsInProfile(), &PtiVGridParProfile, &PtVExtCRSParProfile);
497
498 int const nParProfSize = static_cast<int>(PtiVGridParProfile.size());
499
500 // We have a parallel profile which starts at the coast, but is it long enough to be useful? May have been cut short because it extended outside the grid, or we hit an adjacent profile
501 if (nParProfSize < 3)
502 {
503 // We cannot use this parallel profile, it is too short to calculate along-profile slope, so abandon it and move on to the next parallel profile in this direction
504 nParCoastXLast = nParCoastX;
505 nParCoastYLast = nParCoastY;
506 // LogStream << m_ulIter << ": parallel profile abandoned since too short, starts at [" << nParCoastX << "][" << nParCoastY << "] coastline point " << nThisPointOnCoast << ", length = " << nParProfSize << endl;
507 continue;
508 }
509
510 // This parallel profile is OK, so calculate potential erosion along it. First calculate the length of the parallel profile in external CRS units
511 double const dParProfileLenXY = dGetDistanceBetween(&PtVExtCRSParProfile[0], &PtVExtCRSParProfile[nParProfSize - 1]);
512
513 // Next calculate the distance between profile points, again in external CRS units. Assume that the sample points are equally spaced along the parallel profile (not quite true)
514 double dParSpacingXY = dParProfileLenXY / (nParProfSize - 1);
515
516 // Safety check
517 if (bFPIsEqual(dParSpacingXY, 0.0, TOLERANCE))
518 dParSpacingXY = TOLERANCE;
519
520 // LogStream << "dParSpacingXY = " << dParSpacingXY << endl;
521
522 vector<double> dVParProfileZ(nParProfSize, 0); // Initial (pre-erosion) elevation of both consolidated and unconsolidated sediment for cells 'under' the parallel profile
523 vector<double> dVParProfileDistXY(nParProfSize, 0); // Along-profile distance measured from the coast, in external CRS units
524 vector<double> dVParConsProfileZ(nParProfSize, 0); // Initial (pre-erosion) elevation of consolidated sediment only for cells 'under' the parallel profile
525 vector<double> dVParConsZDiff(nParProfSize, 0);
526 vector<double> dVParConsSlope(nParProfSize, 0);
527
528 for (int i = 0; i < nParProfSize; i++)
529 {
530 int const nXPar = PtiVGridParProfile[i].nGetX();
531 int const nYPar = PtiVGridParProfile[i].nGetY();
532
533 // Is this a sea cell?
534 if (! m_pRasterGrid->m_Cell[nXPar][nYPar].bIsInundated())
535 {
536 // It isn't so move along, nothing to do here
537 // LogStream << m_ulIter << " : [" << nXPar << "][" << nYPar << "] is not inundated" << endl;
538 continue;
539 }
540
541 // Is this cell in a polygon?
542 int const nPolyID = m_pRasterGrid->m_Cell[nXPar][nYPar].nGetPolygonID();
543
544 if (nPolyID == INT_NODATA)
545 {
546 // It isn't. This can happen at the seaward end of polygons TODO 026 Is it a problem?
547 // LogStream << m_ulIter << " : [" << nXPar << "][" << nYPar << "] = {" << dGridCentroidXToExtCRSX(nXPar) << ", " << dGridCentroidYToExtCRSY(nYPar) << "} is not in a polygon" << endl;
548 continue;
549 }
550
551 // Get the number of the highest layer with non-zero thickness
552 int const nTopLayer = m_pRasterGrid->m_Cell[nXPar][nYPar].nGetTopNonZeroLayerAboveBasement();
553
554 // Safety check
555 if (nTopLayer == INT_NODATA)
557
558 if (nTopLayer == NO_NONZERO_THICKNESS_LAYERS)
559 // TODO 025 We are down to basement
560 return RTN_OK;
561
562 // Get the elevation for consolidated sediment only on this cell
563 dVParConsProfileZ[i] = m_pRasterGrid->m_Cell[nXPar][nYPar].dGetConsSedTopForLayerAboveBasement(nTopLayer);
564
565 // Get the elevation for both consolidated and unconsolidated sediment on this cell
566 dVParProfileZ[i] = m_pRasterGrid->m_Cell[nXPar][nYPar].dGetSedimentTopElev();
567
568 // And store the X-Y plane distance from the start of the profile
569 dVParProfileDistXY[i] = i * dParSpacingXY;
570 }
571
572 for (int i = 0; i < nParProfSize - 1; i++)
573 {
574 // For the consolidated-only profile, get the Z differences (already in external CRS units)
575 dVParConsZDiff[i] = dVParConsProfileZ[i] - dVParConsProfileZ[i + 1];
576
577 // Calculate dZ/dXY, the Z slope (i.e. rise over run) in the XY direction. Note that we use the elevation difference on the seaward side of 'this' point
578 dVParConsSlope[i] = dVParConsZDiff[i] / dParSpacingXY;
579 }
580
581 // Sort out the final slope value
582 dVParConsSlope[nParProfSize - 1] = dVParConsSlope[nParProfSize - 2];
583
585 {
586 // Smooth the vector of slopes for the consolidated-only profile
587 dVParConsSlope = dVSmoothProfileSlope(&dVParConsSlope);
588 }
589
590 // Initialize the parallel profile vector with depth / m_dWaveBreakingDepth
591 vector<double> dVParProfileDepthOverDB(nParProfSize, 0); // Depth / wave breaking depth at the parallel profile sample points
592 vector<double> dVParProfileErosionPotential(nParProfSize, 0); // Erosion potential at the parallel profile sample points
593
594 // Calculate the erosion potential along this profile using the shape function which we read in earlier
595 double dTotalErosionPotential = 0;
596
597 // Safety check TODO 061 Is this safety check to depth of breaking a reasonable thing to do?
598 if (dDepthOfBreaking <= 0.0)
599 dDepthOfBreaking = 1e-10;
600
601 for (int i = 0; i < nParProfSize; i++)
602 {
603 // Use the actual depth of water here (i.e. the depth to the top of the unconsolidated sediment, including the thickness of consolidated sediment beneath it)
604 dVParProfileDepthOverDB[i] = m_dThisIterSWL - dVParProfileZ[i];
605 dVParProfileDepthOverDB[i] /= dDepthOfBreaking;
606
607 // Constrain dDepthOverDB[i] to be between 0 (can get small -ve due to rounding errors) and m_dDepthOverDBMax
608 dVParProfileDepthOverDB[i] = tMax(dVParProfileDepthOverDB[i], 0.0);
609 dVParProfileDepthOverDB[i] = tMin(dVParProfileDepthOverDB[i], m_dDepthOverDBMax);
610
611 // And then use the look-up table to find the value of erosion potential at this point on the profile
612 dVParProfileErosionPotential[i] = dLookUpErosionPotential(dVParProfileDepthOverDB[i]);
613
614 // If erosion potential (a -ve value) is tiny, set it to zero
615 if (dVParProfileErosionPotential[i] > -SEDIMENT_ELEV_TOLERANCE)
616 dVParProfileErosionPotential[i] = 0;
617
618 // Keep track of the total erosion potential for this profile
619 dTotalErosionPotential += dVParProfileErosionPotential[i];
620 }
621
622 // Constrain erosion potential at every point on the profile, so that the integral of erosion potential on the whole profile is unity (Walkden and Hall 2005). Note that here, erosion potential is -ve so we must constrain to -1
623 for (int i = 0; i < nParProfSize; i++)
624 {
625 if (dTotalErosionPotential < 0)
626 dVParProfileErosionPotential[i] /= (-dTotalErosionPotential);
627 }
628
629 vector<double> dVParRecessionXY(nParProfSize, 0);
630 vector<double> dVParSCAPEXY(nParProfSize, 0);
631
632 // Calculate recession at every point on the parallel profile
633 for (int i = 0; i < nParProfSize; i++)
634 {
635 // dRecession = dForce * (dBeachProtection / dR) * dErosionPotential * dSlope * dTime
636 // where:
637 // dVRecession [m] is the landward migration distance defined in the profile relative (XY) CRS
638 // dForce is given by Equation 4 in Walkden & Hall, 2005
639 // dVBeachProtection [1] is beach protection factor [1, 0] = [no protection, fully protected] (This is calculated later, see dCalcBeachProtectionFactor())
640 // dVR [m^(9/4)s^(2/3)] is the material strength and some hydrodynamic constant
641 // dVProfileErosionPotential [?] is the erosion potential at each point along the profile
642 // dVSlope [1] is the along-profile slope
643 // m_dTimeStep * 3600 [s] is the time interval in seconds
644 //
645 // dRecession is horizontal recession (along the XY direction):
646 //
647 // dVRecessionXY[i] = (dForce * dVBeachProtection[i] * dVErosionPotentialFunc[i] * dVSlope[i] * m_dTimeStep * 3600) / dVR[i]
648 //
649 //
650 // XY recession must be -ve or zero. If it is +ve then it represents accretion not erosion, which must be described by a different set of equations. So we also need to constrain XY recession to be <= 0
651 dVParRecessionXY[i] = tMin(m_VCoast[nCoast].dGetWaveEnergyAtBreaking(nThisPointOnCoast) * dVParProfileErosionPotential[i] * dVParConsSlope[i] / m_dR, 0.0);
652 dVParSCAPEXY[i] = dVParProfileDistXY[i] - dVParRecessionXY[i];
653
654 // LogStream << m_ulIter << ": [" << nXPar << "][" << nYPar << "] = {" << dGridCentroidXToExtCRSX(nXPar) << ", " << dGridCentroidYToExtCRSY(nYPar) << "} wave energy = " << m_VCoast[nCoast].dGetWaveEnergyAtBreaking(nThisPointOnCoast) << " erosion potential = " << dVParProfileErosionPotential[i] << " slope = " << dVParProfileSlope[i] << " dVParZDiff[i] = " << dVParZDiff[i] << " nParProfSize = " << nParProfSize << endl;
655 }
656
657 vector<double> dVParDeltaZ(nParProfSize, 0);
658
659 // We have calculated the XY-plane recession at every point on the profile, so now convert this to a change in Z-plane elevation at every inundated point on the profile (not the coast point). Again we use the elevation difference on the seaward side of 'this' point
660 for (int i = 1; i < nParProfSize - 1; i++)
661 {
662 // Vertical lowering dZ = dXY * tan(a), where tan(a) is the slope of the SCAPE profile in the XY direction
663 double const dSCAPEHorizDist = dVParSCAPEXY[i + 1] - dVParSCAPEXY[i];
664
665 // Safety check
666 if (bFPIsEqual(dSCAPEHorizDist, 0.0, TOLERANCE))
667 continue;
668
669 double const dSCAPEVertDist = dVParConsProfileZ[i] - dVParConsProfileZ[i + 1];
670 double const dSCAPESlope = dSCAPEVertDist / dSCAPEHorizDist;
671 double dDeltaZ = dVParRecessionXY[i] * dSCAPESlope;
672
673 // Safety check: if thickness model has some jumps, dVConsProfileZ might be very high, limiting dSCAPESlope to 0 because all time erode a high fix quantity
674 if (dSCAPESlope > 1)
675 {
676 dDeltaZ = 0;
677 }
678
679 int const nXPar = PtiVGridParProfile[i].nGetX();
680 int const nYPar = PtiVGridParProfile[i].nGetY();
681
682 // Store the local slope of the consolidated sediment, this is just for output display purposes
683 m_pRasterGrid->m_Cell[nXPar][nYPar].SetLocalConsSlope(dVParConsSlope[i]);
684
685 // dDeltaZ is zero or -ve: if dDeltaZ is zero then do nothing, if -ve then remove some sediment from this cell
686 if (dDeltaZ < 0)
687 {
688 // Has this cell already been eroded during this timestep?
689 if (m_pRasterGrid->m_Cell[nXPar][nYPar].bPotentialPlatformErosion())
690 {
691 // It has
692 double const dPrevPotentialErosion = -m_pRasterGrid->m_Cell[nXPar][nYPar].dGetPotentialPlatformErosion();
693
694 // LogStream << m_ulIter << ": [" << nXPar << "][" << nYPar << "] parallel profile " << nDistFromProfile << " coast points " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast from profile " << nProfile << " has previous potential platform erosion = " << dPrevPotentialErosion << ", current potential platform erosion = " << dDeltaZ << ", max value = " << tMin(dPrevPotentialErosion, dDeltaZ) << endl;
695
696 // Use the larger of the two -ve values
697 dDeltaZ = tMin(dPrevPotentialErosion, dDeltaZ);
698
699 // Adjust this-timestep totals, since this cell has already been eroded
701 m_dThisIterPotentialPlatformErosion += dPrevPotentialErosion; // Since dPrevPotentialErosion is +ve
702 // assert(isfinite(m_dThisIterPotentialPlatformErosion));
703 // assert(m_dThisIterPotentialPlatformErosion >= 0);
704
705 // And also adjust the check values
707 m_dTotPotentialPlatformErosionBetweenProfiles += dPrevPotentialErosion; // Since -ve
708 }
709
710 // Constrain the lowering so we don't get negative slopes or +ve erosion amounts (dDeltaZ must be -ve), this is implicit in SCAPE
711 dDeltaZ = tMax(dDeltaZ, -dVParConsZDiff[i]);
712 dDeltaZ = tMin(dDeltaZ, 0.0);
713 dVParDeltaZ[i] = dDeltaZ;
714
715 // Set the potential (unconstrained) erosion for this cell, it is a +ve value
716 m_pRasterGrid->m_Cell[nXPar][nYPar].SetPotentialPlatformErosion(-dDeltaZ);
717 // LogStream << "[" << nXPar << "][" << nYPar << "] = {" << dGridCentroidXToExtCRSX(nXPar) << ", " << dGridCentroidYToExtCRSY(nYPar) << "} has potential platform erosion = " << -dDeltaZ << endl;
718
719 // Update this-timestep totals
721 m_dThisIterPotentialPlatformErosion -= dDeltaZ; // Since dDeltaZ is a -ve value
722 // assert(isfinite(m_dThisIterPotentialPlatformErosion));
723 // assert(m_dThisIterPotentialPlatformErosion >= 0);
724
725 // Increment the check values
727 m_dTotPotentialPlatformErosionBetweenProfiles -= dDeltaZ; // Since -ve
728 }
729
730 // Finally, calculate the beach protection factor, this will be used in estimating actual (supply-limited) erosion
731 double const dBeachProtectionFactor = dCalcBeachProtectionFactor(nXPar, nYPar, dBreakingWaveHeight);
732 m_pRasterGrid->m_Cell[nXPar][nYPar].SetBeachProtectionFactor(dBeachProtectionFactor);
733 }
734
735 // If desired, save this parallel coastline-normal profile for checking purposes
737 {
738 int const nRet = nSaveParProfile(nCoast, pProfile, nParProfSize, nDirection, nDistFromProfile, &dVParProfileDistXY, &dVParConsProfileZ, &dVParProfileDepthOverDB, &dVParProfileErosionPotential, &dVParConsSlope, &dVParRecessionXY, &dVParDeltaZ, pProfile->pPtiVGetCellsInProfile(), &dVParSCAPEXY);
739
740 if (nRet != RTN_OK)
741 return nRet;
742 }
743
744 // Update for next time round the loop
745 nParCoastXLast = nParCoastX;
746 nParCoastYLast = nParCoastY;
747 }
748
749 return RTN_OK;
750}
751
752//===============================================================================================================================
754//===============================================================================================================================
755void CSimulation::DoActualPlatformErosionOnCell(int const nX, int const nY)
756{
757 // LogStream << m_ulIter << ": doing platform erosion on cell [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
758
759 // Get the beach protection factor, which quantifies the extent to which unconsolidated sediment on the shore platform (beach) protects the shore platform
760 double const dBeachProtectionFactor = m_pRasterGrid->m_Cell[nX][nY].dGetBeachProtectionFactor();
761
762 if (bFPIsEqual(dBeachProtectionFactor, 0.0, TOLERANCE))
763 // The beach is sufficiently thick to prevent any platform erosion (or we are down to basement)
764 return;
765
766 // Get the potential depth of potential erosion, considering beach protection
767 double const dThisPotentialErosion = m_pRasterGrid->m_Cell[nX][nY].dGetPotentialPlatformErosion() * dBeachProtectionFactor;
768
769 // We will be eroding the topmost layer that has non-zero thickness
770 int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
771
772 // Safety check
773 if (nThisLayer == INT_NODATA)
774 {
775 cerr << ERR << "no sediment layer in DoActualPlatformErosionOnCell()" << endl;
776 return;
777 }
778
779 if (nThisLayer == NO_NONZERO_THICKNESS_LAYERS)
780 // No layer with non-zero thickness left, we are down to basement
781 return;
782
783 // OK, we have a layer that can be eroded so find out how much consolidated sediment we have available on this cell
784 double const dExistingAvailableFine = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetFineDepth();
785 double const dExistingAvailableSand = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetSandDepth();
786 double const dExistingAvailableCoarse = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetCoarseDepth();
787
788 // Now partition the total lowering for this cell between the three size fractions: do this by relative erodibility
789 int const nFineWeight = (dExistingAvailableFine > 0 ? 1 : 0);
790 int const nSandWeight = (dExistingAvailableSand > 0 ? 1 : 0);
791 int const nCoarseWeight = (dExistingAvailableCoarse > 0 ? 1 : 0);
792
793 double const dTotErodibility = (nFineWeight * m_dFineErodibilityNormalized) + (nSandWeight * m_dSandErodibilityNormalized) + (nCoarseWeight * m_dCoarseErodibilityNormalized);
794 double dTotActualErosion = 0;
795 double dSandEroded = 0;
796 double dCoarseEroded = 0;
797
798 if (nFineWeight)
799 {
800 // Erode some fine-sized consolidated sediment
801 double const dFineLowering = (m_dFineErodibilityNormalized * dThisPotentialErosion) / dTotErodibility;
802
803 // Make sure we don't get -ve amounts left on the cell
804 double const dFineEroded = tMin(dExistingAvailableFine, dFineLowering);
805 double const dRemaining = dExistingAvailableFine - dFineEroded;
806
807 dTotActualErosion += dFineEroded;
808
809 // Set the value for this layer
810 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetFineDepth(dRemaining);
811
812 // And set the changed-this-timestep switch
813 m_bConsChangedThisIter[nThisLayer] = true;
814
815 // And increment the per-timestep total, also add to the suspended sediment load
818 }
819
820 if (nSandWeight)
821 {
822 // Erode some sand-sized consolidated sediment
823 double const dSandLowering = (m_dSandErodibilityNormalized * dThisPotentialErosion) / dTotErodibility;
824
825 // Make sure we don't get -ve amounts left on the source cell
826 dSandEroded = tMin(dExistingAvailableSand, dSandLowering);
827 double const dRemaining = dExistingAvailableSand - dSandEroded;
828
829 dTotActualErosion += dSandEroded;
830
831 // Set the new value of sand consolidated sediment depth for this layer
832 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetSandDepth(dRemaining);
833
834 // And add this to the depth of sand unconsolidated sediment for this layer
835 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->AddSandDepth(dSandEroded);
836
837 // Set the changed-this-timestep switch
838 m_bConsChangedThisIter[nThisLayer] = true;
839
840 // And increment the per-timestep total
842 }
843
844 if (nCoarseWeight)
845 {
846 // Erode some coarse-sized consolidated sediment
847 double const dCoarseLowering = (m_dCoarseErodibilityNormalized * dThisPotentialErosion) / dTotErodibility;
848
849 // Make sure we don't get -ve amounts left on the source cell
850 dCoarseEroded = tMin(dExistingAvailableCoarse, dCoarseLowering);
851 double const dRemaining = dExistingAvailableCoarse - dCoarseEroded;
852
853 dTotActualErosion += dCoarseEroded;
854
855 // Set the new value of coarse consolidated sediment depth for this layer
856 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->SetCoarseDepth(dRemaining);
857
858 // And add this to the depth of coarse unconsolidated sediment for this layer
859 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->AddCoarseDepth(dCoarseEroded);
860
861 // Set the changed-this-timestep switch
862 m_bConsChangedThisIter[nThisLayer] = true;
863
864 // And increment the per-timestep total
866 }
867
868 // Did we erode anything?
869 if (dTotActualErosion > 0)
870 {
871 // We did, so set the actual erosion value for this cell
872 m_pRasterGrid->m_Cell[nX][nY].SetActualPlatformErosion(dTotActualErosion);
873
874 // Recalculate the elevation of every layer
875 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
876
877 // And update the cell's sea depth
878 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
879
880 // Update per-timestep totals
882
883 // 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)
884 int nPolyID = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
885
886 if (nPolyID == INT_NODATA)
887 {
888 // Can get occasional problems with polygon rasterization near the coastline, so also search the eight adjacent cells
889 array<int, 8> nDirection = {NORTH, NORTH_EAST, EAST, SOUTH_EAST, SOUTH, SOUTH_WEST, WEST, NORTH_WEST};
890 shuffle(nDirection.begin(), nDirection.end(), m_Rand[0]);
891
892 for (int n = 0; n < 8; n++)
893 {
894 int nXAdj;
895 int nYAdj;
896
897 if (nDirection[n] == NORTH)
898 {
899 nXAdj = nX;
900 nYAdj = nY - 1;
901
902 if (nYAdj >= 0)
903 {
904 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
905
906 if (nPolyID != INT_NODATA)
907 break;
908 }
909 }
910
911 else if (nDirection[n] == NORTH_EAST)
912 {
913 nXAdj = nX + 1;
914 nYAdj = nY - 1;
915
916 if ((nXAdj < m_nXGridSize) && (nYAdj >= 0))
917 {
918 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
919
920 if (nPolyID != INT_NODATA)
921 break;
922 }
923 }
924
925 else if (nDirection[n] == EAST)
926 {
927 nXAdj = nX + 1;
928 nYAdj = nY;
929
930 if (nXAdj < m_nXGridSize)
931 {
932 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
933
934 if (nPolyID != INT_NODATA)
935 break;
936 }
937 }
938
939 else if (nDirection[n] == SOUTH_EAST)
940 {
941 nXAdj = nX + 1;
942 nYAdj = nY + 1;
943
944 if ((nXAdj < m_nXGridSize) && (nYAdj < m_nYGridSize))
945 {
946 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
947
948 if (nPolyID != INT_NODATA)
949 break;
950 }
951 }
952
953 else if (nDirection[n] == SOUTH)
954 {
955 nXAdj = nX;
956 nYAdj = nY + 1;
957
958 if (nYAdj < m_nYGridSize)
959 {
960 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
961
962 if (nPolyID != INT_NODATA)
963 break;
964 }
965 }
966
967 else if (nDirection[n] == SOUTH_WEST)
968 {
969 nXAdj = nX - 1;
970 nYAdj = nY + 1;
971
972 if ((nXAdj >= 0) && (nXAdj < m_nXGridSize))
973 {
974 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
975
976 if (nPolyID != INT_NODATA)
977 break;
978 }
979 }
980
981 else if (nDirection[n] == WEST)
982 {
983 nXAdj = nX - 1;
984 nYAdj = nY;
985
986 if (nXAdj >= 0)
987 {
988 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
989
990 if (nPolyID != INT_NODATA)
991 break;
992 }
993 }
994
995 else if (nDirection[n] == NORTH_WEST)
996 {
997 nXAdj = nX - 1;
998 nYAdj = nY - 1;
999
1000 if ((nXAdj >= 0) && (nYAdj >= 0))
1001 {
1002 nPolyID = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetPolygonID();
1003
1004 if (nPolyID != INT_NODATA)
1005 break;
1006 }
1007 }
1008 }
1009 }
1010
1011 // assert(nPolyID < m_VCoast[0].nGetNumPolygons());
1012
1013 // Safety check
1014 if (nPolyID == INT_NODATA)
1015 {
1016 // Uh-oh, we have a problem
1018 LogStream << m_ulIter << ": " << WARN << "platform erosion on cell [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} but this is not in a polygon" << endl;
1019
1020 // m_dDepositionSandDiff and m_dDepositionCoarseDiff are both +ve
1021 m_dDepositionSandDiff += dSandEroded;
1022 m_dDepositionCoarseDiff += dCoarseEroded;
1023
1024 return;
1025 }
1026
1027 // 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)
1028 m_pVCoastPolygon[nPolyID]->AddPlatformErosionUnconsSand(dSandEroded);
1029 m_pVCoastPolygon[nPolyID]->AddPlatformErosionUnconsCoarse(dCoarseEroded);
1030 }
1031}
1032
1033//===============================================================================================================================
1035//===============================================================================================================================
1036bool CSimulation::bCreateErosionPotentialLookUp(vector<double> *VdDepthOverDBIn, vector<double> *VdErosionPotentialIn, vector<double> *VdErosionPotentialFirstDerivIn)
1037{
1038 // Set up a temporary vector to hold the incremental DepthOverDB values
1039 vector<double> VdDepthOverDB;
1040 double dTempDOverDB = 0;
1041
1042 while (dTempDOverDB <= 1.1) // Arbitrary max value, we will adjust this later
1043 {
1044 VdDepthOverDB.push_back(dTempDOverDB); // These are the incremental sample values of DepthOverDB
1045 dTempDOverDB += DEPTH_OVER_DB_INCREMENT;
1046
1047 m_VdErosionPotential.push_back(0); // This will hold the corresponding value of erosion potential for each sample point
1048 }
1049
1050 int const nSize = static_cast<int>(VdDepthOverDB.size());
1051 vector<double> VdDeriv(nSize, 0); // First derivative at the sample points: calculated by the spline function but not subsequently used
1052 vector<double> VdDeriv2(nSize, 0.); // Second derivative at the sample points, ditto
1053 vector<double> VdDeriv3(nSize, 0.); // Third derivative at the sample points, ditto
1054
1055 // Calculate the value of erosion potential (is a -ve value) for each of the sample values of DepthOverDB, and store it for use in the look-up function
1056 hermite_cubic_spline_value(static_cast<int>(VdDepthOverDBIn->size()), &(VdDepthOverDBIn->at(0)), &(VdErosionPotentialIn->at(0)), &(VdErosionPotentialFirstDerivIn->at(0)), nSize, &(VdDepthOverDB[0]), &(m_VdErosionPotential[0]), &(VdDeriv[0]), &(VdDeriv2[0]), &(VdDeriv3[0]));
1057
1058 // Tidy the erosion potential look-up data: cut off values (after the first) for which erosion potential is no longer -ve
1059 int nLastVal = -1;
1060
1061 for (int n = 1; n < nSize - 1; n++)
1062 if (m_VdErosionPotential[n] > 0)
1063 {
1064 nLastVal = n;
1065 break;
1066 }
1067
1068 if (nLastVal > 0)
1069 {
1070 // Erosion potential is no longer -ve at this value of DepthOverDB, so set the maximum value of DepthOverDB that will be used in the simulation (any DepthOverDB value greater than this produces zero erosion potential)
1071 m_dDepthOverDBMax = VdDepthOverDB[nLastVal];
1072 m_VdErosionPotential.erase(m_VdErosionPotential.begin() + nLastVal + 1, m_VdErosionPotential.end());
1073 m_VdErosionPotential.back() = 0;
1074 }
1075
1076 else
1077 // Erosion potential is unbounded, i.e. it is still -ve when we have reached the end of the look-up vector
1078 return false;
1079
1080 // All OK
1081 return true;
1082}
1083
1084//===============================================================================================================================
1086//===============================================================================================================================
1087double CSimulation::dLookUpErosionPotential(double const dDepthOverDB)
1088{
1089 // If dDepthOverDB exceeds the maximum value which we calculated earlier, erosion potential is zero
1090 if (dDepthOverDB > m_dDepthOverDBMax)
1091 return 0;
1092
1093 // OK, dDepthOverDB is less than the maximum so look up a corresponding value for erosion potential. The look-up index is dDepthOverDB divided by (the Depth Over DB increment used when creating the look-up vector). But since this look-up index may not be an integer, split the look-up index into integer and fractional parts and deal with each separately
1094 double const dErosionPotential = dGetInterpolatedValue(&m_VdDepthOverDB, &m_VdErosionPotential, dDepthOverDB, false);
1095
1096 return dErosionPotential;
1097}
1098
1099//===============================================================================================================================
1101//===============================================================================================================================
1102double CSimulation::dCalcBeachProtectionFactor(int const nX, int const nY, double const dBreakingWaveHeight)
1103{
1104 // Safety check
1105 if (bFPIsEqual(dBreakingWaveHeight, DBL_NODATA, TOLERANCE))
1106 return 0;
1107
1108 // We are considering the unconsolidated sediment (beach) of the topmost layer that has non-zero thickness
1109 int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1110
1111 // Safety check
1112 if (nThisLayer == INT_NODATA)
1113 {
1114 cerr << ERR << "no sediment layer in dCalcBeachProtectionFactor()" << endl;
1115 return 0;
1116 }
1117
1118 if (nThisLayer == NO_NONZERO_THICKNESS_LAYERS)
1119 // There are no layers with non-zero thickness left (i.e. we are down to basement) so no beach protection
1120 return 0;
1121
1122 // In SCAPE, 0.23 * the significant breaking wave height is assumed to be the maximum depth of beach that waves can penetrate to erode a platform. For depths less than this, the beach protective ability is assumed to vary linearly
1123 double const dBeachDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->dGetUnconsolidatedThickness();
1124 double const dMaxPenetrationDepth = BEACH_PROTECTION_HB_RATIO * dBreakingWaveHeight;
1125 double dFactor = 0;
1126
1127 if (dMaxPenetrationDepth > 0)
1128 dFactor = tMax(1 - (dBeachDepth / dMaxPenetrationDepth), 0.0);
1129
1130 // LogStream << m_ulIter << ": cell[" << nX << "][" << nY << "] has beach protection factor = " << dFactor << endl;
1131
1132 return dFactor;
1133}
1134
1135//===============================================================================================================================
1137//===============================================================================================================================
1139{
1140 for (int nX = 0; nX < m_nXGridSize; nX++)
1141 {
1142 for (int nY = 0; nY < m_nYGridSize; nY++)
1143 {
1144 if ((m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea()) && (bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetBeachProtectionFactor(), DBL_NODATA, TOLERANCE)))
1145 {
1146 // This is a sea cell, and it has an initialized beach protection value. So look at its N-S and W-E neighbours
1147 int nXTmp;
1148 int nYTmp;
1149 int nAdjacent = 0;
1150 double dBeachProtection = 0;
1151
1152 // North
1153 nXTmp = nX;
1154 nYTmp = nY - 1;
1155
1156 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor(), DBL_NODATA, TOLERANCE)))
1157 {
1158 nAdjacent++;
1159 dBeachProtection += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor();
1160 }
1161
1162 // East
1163 nXTmp = nX + 1;
1164 nYTmp = nY;
1165
1166 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor(), DBL_NODATA, TOLERANCE)))
1167 {
1168 nAdjacent++;
1169 dBeachProtection += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor();
1170 }
1171
1172 // South
1173 nXTmp = nX;
1174 nYTmp = nY + 1;
1175
1176 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor(), DBL_NODATA, TOLERANCE)))
1177 {
1178 nAdjacent++;
1179 dBeachProtection += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor();
1180 }
1181
1182 // West
1183 nXTmp = nX - 1;
1184 nYTmp = nY;
1185
1186 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor(), DBL_NODATA, TOLERANCE)))
1187 {
1188 nAdjacent++;
1189 dBeachProtection += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetBeachProtectionFactor();
1190 }
1191
1192 // If this sea cell has four neighbours with initialized beach protection values, then assume that it should not have an uninitialized beach protection value. Set it to the average of its neighbours
1193 if (nAdjacent == 4)
1194 {
1195 m_pRasterGrid->m_Cell[nX][nY].SetBeachProtectionFactor(dBeachProtection / 4);
1196 }
1197 }
1198 }
1199 }
1200}
1201
1202//===============================================================================================================================
1204//===============================================================================================================================
1206{
1207 for (int nX = 0; nX < m_nXGridSize; nX++)
1208 {
1209 for (int nY = 0; nY < m_nYGridSize; nY++)
1210 {
1211 if ((m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea()) && (bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetPotentialPlatformErosion(), 0.0, TOLERANCE)))
1212 {
1213 // This is a sea cell, it has a zero potential platform erosion value. So look at its N-S and W-E neighbours
1214 int nXTmp;
1215 int nYTmp;
1216 int nAdjacent = 0;
1217 double dPotentialPlatformErosion = 0;
1218
1219 // North
1220 nXTmp = nX;
1221 nYTmp = nY - 1;
1222
1223 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion(), 0.0, TOLERANCE)))
1224 {
1225 nAdjacent++;
1226 dPotentialPlatformErosion += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion();
1227 }
1228
1229 // East
1230 nXTmp = nX + 1;
1231 nYTmp = nY;
1232
1233 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion(), 0.0, TOLERANCE)))
1234 {
1235 nAdjacent++;
1236 dPotentialPlatformErosion += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion();
1237 }
1238
1239 // South
1240 nXTmp = nX;
1241 nYTmp = nY + 1;
1242
1243 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion(), 0.0, TOLERANCE)))
1244 {
1245 nAdjacent++;
1246 dPotentialPlatformErosion += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion();
1247 }
1248
1249 // West
1250 nXTmp = nX - 1;
1251 nYTmp = nY;
1252
1253 if ((bIsWithinValidGrid(nXTmp, nYTmp)) && (! bFPIsEqual(m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion(), 0.0, TOLERANCE)))
1254 {
1255 nAdjacent++;
1256 dPotentialPlatformErosion += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetPotentialPlatformErosion();
1257 }
1258
1259 // If this sea cell has four neighbours with non-zero potential platform erosion values, then assume that it should not have a zero potential platform erosion value. Set it to the average of its neighbours
1260 if (nAdjacent == 4)
1261 {
1262 double const dThisPotentialPlatformErosion = dPotentialPlatformErosion / 4;
1263
1264 m_pRasterGrid->m_Cell[nX][nY].SetPotentialPlatformErosion(dThisPotentialPlatformErosion);
1265
1266 // Update this-timestep totals
1268 m_dThisIterPotentialPlatformErosion += dThisPotentialPlatformErosion;
1269 // assert(isfinite(m_dThisIterPotentialPlatformErosion));
1270
1271 // Increment the check values
1273 m_dTotPotentialPlatformErosionBetweenProfiles += dThisPotentialPlatformErosion;
1274 }
1275 }
1276 }
1277 }
1278}
1279
1280//===============================================================================================================================
1282//===============================================================================================================================
1283void CSimulation::ConstructParallelProfile(int const nProfileStartX, int const nProfileStartY, int const nParCoastX, int const nParCoastY, int const nProfSize, vector<CGeom2DIPoint> *const pPtViGridProfile, vector<CGeom2DIPoint> *pPtiVGridParProfile, vector<CGeom2DPoint> *pPtVExtCRSParProfile)
1284{
1285 // OK, we have the coastline start point for the parallel profile. Now construct a temporary profile, parallel to the coastline-normal profile, starting from this point
1286 int const nXOffset = nParCoastX - nProfileStartX;
1287 int const nYOffset = nParCoastY - nProfileStartY;
1288
1289 // Append co-ord values for the temporary profile
1290 for (int nProfileStartPoint = 0; nProfileStartPoint < nProfSize; nProfileStartPoint++)
1291 {
1292 // Get the grid coordinates of the cell which is this distance seaward, from the coastline-normal profile
1293 int const nXProf = pPtViGridProfile->at(nProfileStartPoint).nGetX();
1294 int const nYProf = pPtViGridProfile->at(nProfileStartPoint).nGetY();
1295
1296 // Now calculate the grid coordinates of this cell, which is potentially in the parallel profile
1297 int const nXPar = nXProf + nXOffset;
1298 int const nYPar = nYProf + nYOffset;
1299
1300 // Is this cell within the grid? If not, cut short the profile
1301 if (! bIsWithinValidGrid(nXPar, nYPar))
1302 {
1303 // LogStream << "NOT WITHIN GRID [" << nXPar << "][" << nYPar << "]" << endl;
1304 return;
1305 }
1306
1307 // Have we hit an adjacent coastline-normal profile? If so, cut short
1308 if (m_pRasterGrid->m_Cell[nXPar][nYPar].bIsProfile())
1309 {
1310 // LogStream << "HIT PROFILE " << m_pRasterGrid->m_Cell[nXPar][nYPar].nGetProfileID() << " at [" << nXPar << "][" << nYPar << "] = {" << dGridCentroidXToExtCRSX(nXPar) << ", " << dGridCentroidYToExtCRSY(nYPar) << "}" << endl;
1311 return;
1312 }
1313
1314 // OK, append the cell details
1315 pPtiVGridParProfile->push_back(CGeom2DIPoint(nXPar, nYPar));
1316 pPtVExtCRSParProfile->push_back(CGeom2DPoint(dGridCentroidXToExtCRSX(nXPar), dGridCentroidYToExtCRSY(nYPar)));
1317 }
1318}
Contains CGeom2DIPoint definitions.
Geometry class used to represent 2D point objects with integer coordinates.
Definition 2di_point.h:29
Geometry class used to represent 2D point objects with floating-point coordinates.
Definition 2d_point.h:27
double dGetY(void) const
Returns the CGeom2DPoint object's double Y coordinate.
Definition 2d_point.cpp:53
double dGetX(void) const
Returns the CGeom2DPoint object's double X coordinate.
Definition 2d_point.cpp:47
Geometry class used to represent coast profile objects.
Definition profile.h:37
int nGetProfileID(void) const
Returns the profile's this-coast ID.
Definition profile.cpp:80
int nGetCoastPoint(void) const
Returns the coast point at which the profile starts.
Definition profile.cpp:86
int nGetNumCellsInProfile(void) const
Returns the number of cells in the profile.
Definition profile.cpp:535
CGeom2DPoint * pPtGetPointInProfile(int const)
Returns a single point (external CRS) from the profile.
Definition profile.cpp:367
bool bOKIncStartAndEndOfCoast(void) const
Returns true if this is a problem-free profile (however it could be a start-of-coast or an end-of-coa...
Definition profile.cpp:265
int nGetProfileSize(void) const
Returns the number of external CRS points in the profile (only two, initally; and always just two for...
Definition profile.cpp:361
vector< CGeom2DIPoint > * pPtiVGetCellsInProfile(void)
Returns all cells (grid CRS) in the profile.
Definition profile.cpp:515
int m_nLogFileDetail
The level of detail in the log file output. Can be LOG_FILE_LOW_DETAIL, LOG_FILE_MIDDLE_DETAIL,...
Definition simulation.h:574
double m_dThisIterSWL
The still water level for this timestep (this includes tidal changes and any long-term SWL change)
Definition simulation.h:712
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
int m_nXGridSize
The size of the grid in the x direction.
Definition simulation.h:457
ofstream LogStream
vector< CRWCoast > m_VCoast
The coastline objects.
double m_dFineErodibilityNormalized
Relative erodibility of fine unconsolidated beach sediment, normalized.
Definition simulation.h:796
double m_dThisIterActualPlatformErosionCoarseCons
Total actual platform erosion (coarse consolidated sediment) for this iteration (depth in m)
Definition simulation.h:841
int nDoAllShorePlatFormErosion(void)
Does platform erosion on all coastlines by first calculating platform erosion on coastline-normal pro...
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:460
vector< double > m_VdErosionPotential
For erosion potential lookup.
void ConstructParallelProfile(int const, int const, int const, int const, int const, vector< CGeom2DIPoint > *const, vector< CGeom2DIPoint > *, vector< CGeom2DPoint > *)
Constructs a parallel coastline-normal profile.
int nSaveParProfile(int const, CGeomProfile const *, int const, int const, int const, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< CGeom2DIPoint > *const, vector< double > const *) const
Save a coastline-normal parallel profile.
static double dGetDistanceBetween(CGeom2DPoint const *, CGeom2DPoint const *)
Returns the distance (in external CRS) between two points.
double m_dCoarseErodibilityNormalized
Relative erodibility of coarse unconsolidated beach sediment, normalized.
Definition simulation.h:802
vector< double > m_VdDepthOverDB
For erosion potential lookup.
void FillInBeachProtectionHoles(void)
Fills in 'holes' in the beach protection i.e. orphan cells which get omitted because of rounding prob...
double m_dTotPotentialPlatformErosionBetweenProfiles
Total potential platform erosion between profiles.
Definition simulation.h:898
double m_dDepositionCoarseDiff
Error term: if we are unable to deposit enough unconslidated coarse on polygon(s),...
Definition simulation.h:889
unsigned long m_ulThisIterNumActualPlatformErosionCells
The number of grid cells on which actual platform erosion occurs, for this iteration.
Definition simulation.h:619
int nCalcPotentialPlatformErosionBetweenProfiles(int const, CGeomProfile *, int const)
Calculates potential platform erosion on cells to one side of a given coastline-normal profile,...
double m_dDepositionSandDiff
Error term: if we are unable to deposit enough unconslidated sand on polygon(s), this is held over to...
Definition simulation.h:886
int nSaveProfile(int const, CGeomProfile const *, int const, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *, vector< CGeom2DIPoint > *const, vector< double > const *) const
Save a coastline-normal profile.
double m_dThisIterPotentialPlatformErosion
Total potential platform erosion (all size classes of consolidated sediment) for this iteration (dept...
Definition simulation.h:832
vector< double > dVSmoothProfileSlope(vector< double > *) const
Does running-mean smoothing of the slope of a coastline-normal profile.
double m_dDepthOverDBMax
Maximum value of deoth over DB, is used in erosion potential look-up function.
Definition simulation.h:892
double m_dSandErodibilityNormalized
Relative erodibility of sand unconsolidated beach sediment, normalized.
Definition simulation.h:799
double m_dR
Coast platform resistance to erosion R, see Walkden & Hall, 2011.
Definition simulation.h:769
double m_dThisIterActualPlatformErosionSandCons
Total actual platform erosion (sand consolidated sediment) for this iteration (depth in m)
Definition simulation.h:838
double dGridCentroidYToExtCRSY(int const) const
Given the integer Y-axis ordinate of a cell in the raster grid CRS, returns the external CRS Y-axis o...
Definition gis_utils.cpp:78
double m_dTotPotentialPlatformErosionOnProfiles
Total potential platform erosion on profiles.
Definition simulation.h:895
int nCalcPotentialPlatformErosionOnProfile(int const, CGeomProfile *)
Calculates potential (i.e. unconstrained by available sediment) erosional lowering of the shore platf...
void DoActualPlatformErosionOnCell(int const, int const)
Calculates actual (constrained by available sediment) erosion of the consolidated shore platform on a...
double m_dThisIterFineSedimentToSuspension
Total fine unconsolidated sediment in suspension for this iteration (depth in m)
Definition simulation.h:862
void FillPotentialPlatformErosionHoles(void)
Fills in 'holes' in the potential platform erosion i.e. orphan cells which get omitted because of rou...
double dLookUpErosionPotential(double const)
The erosion potential lookup: it returns a value for erosion potential given a value of Depth Over DB...
unsigned long m_ulThisIterNumPotentialPlatformErosionCells
The number of grid cells on which potential platform erosion occurs, for this iteration.
Definition simulation.h:616
default_random_engine m_Rand[NUMBER_OF_RNGS]
The c++11 random number generators.
unsigned long m_ulTotPotentialPlatformErosionOnProfiles
The number of cells on which on-profile average potential shore platform erosion occurs.
Definition simulation.h:631
double m_dCellArea
Area of a cell (in external CRS units)
Definition simulation.h:661
unsigned long m_ulTotPotentialPlatformErosionBetweenProfiles
The number of cells on which between-profile average potential shore platform erosion occurs.
Definition simulation.h:634
bool bIsWithinValidGrid(int const, int const) const
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:598
double dGridCentroidXToExtCRSX(int const) const
Definition gis_utils.cpp:68
double dGetInterpolatedValue(vector< double > const *, vector< double > const *, double, bool)
double dCalcBeachProtectionFactor(int const, int const, double const)
Calculates the (inverse) beach protection factor as in SCAPE: 0 is fully protected,...
bool m_bOutputProfileData
Output profile data?
Definition simulation.h:337
bool bCreateErosionPotentialLookUp(vector< double > *, vector< double > *, vector< double > *)
Creates a look-up table for erosion potential, given depth over DB.
double m_dThisIterActualPlatformErosionFineCons
Total actual platform erosion (fine consolidated sediment) for this iteration (depth in m)
Definition simulation.h:835
int m_nProfileSmoothWindow
The size of the window used for running-mean coast-normal profile smoothing (must be odd)
Definition simulation.h:487
vector< bool > m_bConsChangedThisIter
One element per layer: has the consolidated sediment of this layer been changed during this iteration...
bool m_bOutputParallelProfileData
Output parallel profile data?
Definition simulation.h:340
vector< CGeomCoastPolygon * > m_pVCoastPolygon
Pointers to coast polygon objects, in down-coast sequence TODO 044 Will need to use global polygon ID...
This file contains global definitions for CoastalME.
int const NO_NONZERO_THICKNESS_LAYERS
Definition cme.h:778
int const INT_NODATA
Definition cme.h:476
T tMin(T a, T b)
Definition cme.h:1265
double const TOLERANCE
Definition cme.h:823
int const RTN_ERR_NO_TOP_LAYER
Definition cme.h:738
string const ERR
Definition cme.h:902
int const DIRECTION_DOWNCOAST
Definition cme.h:506
int const SOUTH_EAST
Definition cme.h:500
double const DEPTH_OVER_DB_INCREMENT
Definition cme.h:815
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:491
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1303
int const SOUTH_WEST
Definition cme.h:502
T tMax(T a, T b)
Definition cme.h:1252
int const NORTH_WEST
Definition cme.h:504
int const LOG_FILE_HIGH_DETAIL
Definition cme.h:492
string const WARN
Definition cme.h:903
double const SEDIMENT_ELEV_TOLERANCE
Definition cme.h:825
int const SOUTH
Definition cme.h:501
int const LOG_FILE_ALL
Definition cme.h:493
int const EAST
Definition cme.h:499
int const RTN_OK
Definition cme.h:694
int const NORTH_EAST
Definition cme.h:498
int const DIRECTION_UPCOAST
Definition cme.h:507
double const DBL_NODATA
Definition cme.h:834
double const BEACH_PROTECTION_HB_RATIO
Definition cme.h:811
int const NORTH
Definition cme.h:497
int const WEST
Definition cme.h:503
Contains CRWCoast definitions.
void hermite_cubic_spline_value(int const nn, double *const xn, double *const fn, double *const dn, int const n, double *const x, double *f, double *d, double *s, double *t)
This is part of a C++ library from http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_cubic/hermite_...
Definitions of some routines from the hermite_cubic library.
Contains CSimulation definitions.