CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
init_grid.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 <cstdio>
27#include <climits>
28
29#include <iostream>
30using std::cerr;
31using std::endl;
32
33#ifdef _OPENMP
34#include <omp.h>
35#endif
36
37#include "cme.h"
38#include "cell.h"
39#include "coast.h"
40#include "simulation.h"
41#include "raster_grid.h"
42
43//===============================================================================================================================
45//===============================================================================================================================
47{
48 // Clear all vector coastlines, profiles, and polygons
49 for (int i = 0; i < static_cast<int>(m_pVCoastPolygon.size()); i++)
50 delete m_pVCoastPolygon[i];
51
52 m_pVCoastPolygon.clear();
53 m_VCoast.clear();
54
55 // m_VFloodWaveSetup.clear();
58
59 // Do some every-timestep initialization
60 m_nXMinBoundingBox = INT_MAX;
61 m_nXMaxBoundingBox = INT_MIN;
62 m_nYMinBoundingBox = INT_MAX;
63 m_nYMaxBoundingBox = INT_MIN;
64
69
73
95 m_dThisIterLeftGridUnconsFine = // TODO 067 Suspended fine sediment never decreases i.e. no suspended fine sediment ever leaves the grid. Is this OK?
101
102 for (int n = 0; n < m_nLayers; n++)
103 {
104 m_bConsChangedThisIter[n] = false;
105 m_bUnconsChangedThisIter[n] = false;
106 }
107
108 // Re-calculate the depth of closure, in case deep water wave properties have changed
110
111 int nZeroThickness = 0;
112
121
122 // And go through all cells in the RasterGrid array
123 // Use OpenMP parallel loop with reduction clauses for thread-safe accumulation
124#ifdef _OPENMP
125#pragma omp parallel for collapse(2) reduction(+ : nZeroThickness) \
126 reduction(+ : m_dStartIterConsFineAllCells, m_dStartIterConsSandAllCells, m_dStartIterConsCoarseAllCells) \
127 reduction(+ : m_dStartIterSuspFineAllCells, m_dStartIterUnconsFineAllCells, m_dStartIterUnconsSandAllCells, m_dStartIterUnconsCoarseAllCells)
128#endif
129
130 for (int nX = 0; nX < m_nXGridSize; nX++)
131 {
132 for (int nY = 0; nY < m_nYGridSize; nY++)
133 {
134 // Re-initialize values for this cell
135 m_pRasterGrid->m_Cell[nX][nY].InitCell();
136
137 if (m_ulIter == 1)
138 {
139 // For the first timestep only, check to see that all cells have some sediment on them
140 double const dSedThickness = m_pRasterGrid->m_Cell[nX][nY].dGetTotAllSedThickness();
141
142 if (dSedThickness <= 0)
143 {
144 nZeroThickness++;
145
146 // Note: Logging from parallel regions can cause race conditions, but this is for debugging only
147 // In production, consider collecting problematic cells and logging after the parallel region
149 {
150#ifdef _OPENMP
151#pragma omp critical(logging)
152#endif
153 LogStream << m_ulIter << ": " << WARN << "total sediment thickness is " << dSedThickness << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
154 }
155 }
156
157 // For the first timestep only, calculate the elevation of all this cell's layers. During the rest of the simulation, each cell's elevation is re-calculated just after any change occurs on that cell
158 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
159 }
160
161 // Note that these totals include sediment which is both within and outside the polygons (because we have not yet defined polygons for this iteration, duh!)
162 m_dStartIterConsFineAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetTotConsFineThickConsiderNotch();
163 m_dStartIterConsSandAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetTotConsSandThickConsiderNotch();
164 m_dStartIterConsCoarseAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetTotConsCoarseThickConsiderNotch();
165
166 m_dStartIterSuspFineAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetSuspendedSediment();
167 m_dStartIterUnconsFineAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetTotUnconsFine();
168 m_dStartIterUnconsSandAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetTotUnconsSand();
169 m_dStartIterUnconsCoarseAllCells += m_pRasterGrid->m_Cell[nX][nY].dGetTotUnconsCoarse();
170
172 {
173 // If we have just a single measurement for deep water waves (either given by the user, or from a single wave station) then set all cells, even dry land cells, to the same value for deep water wave height, deep water wave orientation, and deep water period
174 m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWaveHeight(m_dAllCellsDeepWaterWaveHeight);
175 m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWaveAngle(m_dAllCellsDeepWaterWaveAngle);
176 m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWavePeriod(m_dAllCellsDeepWaterWavePeriod);
177 }
178 }
179 }
180
182 {
183 // Each cell's value for deep water wave height and deep water wave orientation is interpolated from multiple user-supplied values
184 int const nRet = nInterpolateAllDeepWaterWaveValues();
185
186 if (nRet != RTN_OK)
187 return nRet;
188
189 /* for (int n = 0; n < m_VlDeepWaterWaveValuesAtTimestep.size(); n++)
190 {
191 if (m_ulIter == m_VlDeepWaterWaveValuesAtTimestep[n])
192 {
193 // OK, this timestep we are doing the calculation
194 if (m_VlDeepWaterWaveValuesAtTimestep[n] > 1)
195 {
196 // TODO 036 For every timestep after the first, read in new values before doing the interpolation
197 }
198
199 // Interpolate values each cell's values for deep water height and orientation from user-supplied values
200 int nRet = nInterpolateAllDeepWaterWaveValues();
201 if (nRet != RTN_OK)
202 return nRet;
203 }
204 }*/
205 }
206
207 if (nZeroThickness > 0)
208 {
209 cerr << m_ulIter << ": " << WARN << nZeroThickness << (nZeroThickness > 1 ? " cells" : " cell") << "( out of " << m_nXGridSize * m_nYGridSize << ") have no sediment, is this correct?" << endl;
210 LogStream << m_ulIter << ": " << WARN << nZeroThickness << " cells have no sediment, is this correct?" << endl;
211 }
212
213 // // DEBUG CODE ===========================================================================================================
214 // string strOutFile = m_strOutPath;
215 // strOutFile += "init_deep_water_wave_height_";
216 // strOutFile += to_string(m_ulIter);
217 // strOutFile += ".tif";
218 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
219 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
220 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
221 // pDataSet->SetGeoTransform(m_dGeoTransform);
222 // double* pdRaster = new double[m_ulNumCells];
223 // int nn = 0;
224 // for (int nY = 0; nY < m_nYGridSize; nY++)
225 // {
226 // for (int nX = 0; nX < m_nXGridSize; nX++)
227 // {
228 // // Write this value to the array
229 // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveHeight();
230 // nn++;
231 // }
232 // }
233 //
234 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
235 // pBand->SetNoDataValue(m_nMissingValue);
236 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
237 //
238 // if (nRet == CE_Failure)
239 // return RTN_ERR_GRIDCREATE;
240 //
241 // GDALClose(pDataSet);
242 // // DEBUG CODE ===========================================================================================================
243
244 // // DEBUG CODE ===========================================================================================================
245 // strOutFile = m_strOutPath;
246 // strOutFile += "init_deep_water_wave_angle_";
247 // strOutFile += to_string(m_ulIter);
248 // strOutFile += ".tif";
249 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
250 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
251 // pDataSet->SetGeoTransform(m_dGeoTransform);
252 // nn = 0;
253 // for (int nY = 0; nY < m_nYGridSize; nY++)
254 // {
255 // for (int nX = 0; nX < m_nXGridSize; nX++)
256 // {
257 // // Write this value to the array
258 // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveAngle();
259 // nn++;
260 // }
261 // }
262 //
263 // pBand = pDataSet->GetRasterBand(1);
264 // pBand->SetNoDataValue(m_nMissingValue);
265 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
266 //
267 // if (nRet == CE_Failure)
268 // return RTN_ERR_GRIDCREATE;
269 //
270 // GDALClose(pDataSet);
271 // // DEBUG CODE ===========================================================================================================
272
273 // // DEBUG CODE ===========================================================================================================
274 // strOutFile = m_strOutPath;
275 // strOutFile += "init_water_wave_angle_";
276 // strOutFile += to_string(m_ulIter);
277 // strOutFile += ".tif";
278 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
279 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
280 // pDataSet->SetGeoTransform(m_dGeoTransform);
281 // nn = 0;
282 // for (int nY = 0; nY < m_nYGridSize; nY++)
283 // {
284 // for (int nX = 0; nX < m_nXGridSize; nX++)
285 // {
286 // // Write this value to the array
287 // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
288 // nn++;
289 // }
290 // }
291 //
292 // pBand = pDataSet->GetRasterBand(1);
293 // pBand->SetNoDataValue(m_nMissingValue);
294 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
295 //
296 // if (nRet == CE_Failure)
297 // return RTN_ERR_GRIDCREATE;
298 //
299 // GDALClose(pDataSet);
300 // // DEBUG CODE ===========================================================================================================
301
302 // // DEBUG CODE ===========================================================================================================
303 // strOutFile = m_strOutPath;
304 // strOutFile += "init_water_wave_height_";
305 // strOutFile += to_string(m_ulIter);
306 // strOutFile += ".tif";
307 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
308 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
309 // pDataSet->SetGeoTransform(m_dGeoTransform);
310 // nn = 0;
311 // for (int nY = 0; nY < m_nYGridSize; nY++)
312 // {
313 // for (int nX = 0; nX < m_nXGridSize; nX++)
314 // {
315 // // Write this value to the array
316 // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
317 // nn++;
318 // }
319 // }
320 //
321 // pBand = pDataSet->GetRasterBand(1);
322 // pBand->SetNoDataValue(m_nMissingValue);
323 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
324 //
325 // if (nRet == CE_Failure)
326 // return RTN_ERR_GRIDCREATE;
327 //
328 // GDALClose(pDataSet);
329 // delete[] pdRaster;
330 // // DEBUG CODE ===========================================================================================================
331
332 return RTN_OK;
333}
Contains CGeomCell definitions.
double m_dThisIterPotentialBeachErosion
Total potential beach erosion (all size classes of unconsolidated sediment) for this iteration (depth...
Definition simulation.h:844
void CalcDepthOfClosure(void)
Calculate the depth of closure.
Definition utils.cpp:2648
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_dAllCellsDeepWaterWaveHeight
Deep water wave height (m) for all sea cells.
Definition simulation.h:754
int m_nYMinBoundingBox
The minimum y value of the bounding box.
Definition simulation.h:544
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
double m_dThisIterBeachErosionCoarse
Total actual beach erosion (coarse unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:853
vector< CRWCoast > m_VFloodWaveSetupSurgeRunup
TODO 007 Info needed.
double m_dStartIterUnconsCoarseAllCells
Depth (m) of coarse unconsolidated sediment at the start of the simulation, all cells (both inside an...
Definition simulation.h:994
vector< CRWCoast > m_VCoast
The coastline objects.
bool m_bSingleDeepWaterWaveValues
Do we have just a point source for (i.e. only a single measurement of) deep water wave values.
Definition simulation.h:385
double m_dThisIterActualPlatformErosionCoarseCons
Total actual platform erosion (coarse consolidated sediment) for this iteration (depth in m)
Definition simulation.h:841
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:460
double m_dStartIterSuspFineInPolygons
Depth (m) of fine suspended sediment at the start of the simulation (only cells in polygons)
Definition simulation.h:985
unsigned long m_ulThisIterNumCoastCells
The number of grid cells which are marked as coast, for this iteration.
Definition simulation.h:613
double m_dThisIterCliffCollapseErosionCoarseUncons
This-iteration total of coarse unconsolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:934
vector< CRWCoast > m_VFloodWaveSetupSurge
TODO 007 Info needed.
unsigned long m_ulThisIterNumActualPlatformErosionCells
The number of grid cells on which actual platform erosion occurs, for this iteration.
Definition simulation.h:619
int m_nXMaxBoundingBox
The maximum x value of the bounding box.
Definition simulation.h:541
double m_dThisIterCliffCollapseErosionFineUncons
This-iteration total of fine unconsolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:928
int m_nLayers
The number of sediment layers.
Definition simulation.h:463
double m_dStartIterConsCoarseAllCells
Depth (m) of coarse consolidated sediment at the start of the simulation, all cells (both inside and ...
double m_dStartIterConsSandAllCells
Depth (m) of sand consolidated sediment at the start of the simulation, all cells (both inside and ou...
unsigned long m_ulThisIterNumBeachDepositionCells
The number of grid cells on which beach (unconsolidated sediment) deposition occurs,...
Definition simulation.h:628
double m_dThisiterUnconsCoarseInput
Depth (m) of coarse unconsolidated sediment added, at this iteration.
Definition simulation.h:979
double m_dThisIterPotentialPlatformErosion
Total potential platform erosion (all size classes of consolidated sediment) for this iteration (dept...
Definition simulation.h:832
double m_dThisiterUnconsFineInput
Depth (m) of fine unconsolidated sediment added, at this iteration.
Definition simulation.h:973
bool m_bHaveWaveStationData
Do we have wave station data?
Definition simulation.h:388
double m_dThisIterActualPlatformErosionSandCons
Total actual platform erosion (sand consolidated sediment) for this iteration (depth in m)
Definition simulation.h:838
double m_dThisIterBeachDepositionCoarse
Total beach deposition (coarse unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:859
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_dAllCellsDeepWaterWaveAngle
Deep water wave angle for all sea cells.
Definition simulation.h:757
double m_dThisIterUnconsCoarseCliffDeposition
This-iteration total of coarse unconsolidated sediment deposited due to cliff collapse (m^3)
Definition simulation.h:949
double m_dThisIterTotSeaDepth
Total sea depth (m) for this iteration.
Definition simulation.h:829
double m_dThisIterUnconsSandCliffDeposition
This-iteration total of sand unconsolidated sediment deposited due to cliff collapse (m^3)
Definition simulation.h:946
int nInterpolateAllDeepWaterWaveValues(void)
If the user supplies multiple deep water wave height and angle values, this routine interplates these...
double m_dThisIterPotentialSedLostBeachErosion
Total unconsolidated sediment from beach erosion (all size classes) lost from the grid this iteration...
Definition simulation.h:865
double m_dStartIterConsFineAllCells
Depth (m) of fine consolidated sediment at the start of the simulation, all cells (both inside and ou...
Definition simulation.h:997
double m_dThisIterBeachErosionSand
Total actual beach erosion (sand unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:850
double m_dThisIterFineSedimentToSuspension
Total fine unconsolidated sediment in suspension for this iteration (depth in m)
Definition simulation.h:862
double m_dThisIterBeachErosionFine
Total actual beach erosion (fine unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:847
unsigned long m_ulThisIterNumActualBeachErosionCells
The number of grid cells on which actual beach (unconsolidated sediment) erosion occurs,...
Definition simulation.h:625
int m_nYMaxBoundingBox
The maximum y value of the bounding box.
Definition simulation.h:547
vector< bool > m_bUnconsChangedThisIter
One element per layer: has the consolidated sediment of this layer been changed during this iteration...
unsigned long m_ulThisIterNumPotentialPlatformErosionCells
The number of grid cells on which potential platform erosion occurs, for this iteration.
Definition simulation.h:616
double m_dStartIterUnconsSandAllCells
Depth (m) of sand unconsolidated sediment at the start of the simulation, all cells (both inside and ...
Definition simulation.h:991
double m_dThisIterCliffCollapseErosionFineCons
This-iteration total of fine consolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:937
unsigned long m_ulThisIterNumSeaCells
The number of grid cells which are marked as sea, for this iteration.
Definition simulation.h:610
double m_dStartIterUnconsFineAllCells
Depth (m) of fine unconsolidated sediment at the start of the simulation, all cells (both inside and ...
Definition simulation.h:988
int m_nXMinBoundingBox
The minimum x value of the bounding box.
Definition simulation.h:538
unsigned long m_ulThisIterNumPotentialBeachErosionCells
The number of grid cells on which potential beach (unconsolidated sediment) erosion occurs,...
Definition simulation.h:622
double m_dThisIterCliffCollapseErosionSandCons
This-iteration total of sand consolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:940
double m_dThisIterBeachDepositionSand
Total beach deposition (sand unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:856
double m_dThisIterLeftGridUnconsCoarse
Total coarse unconsolidated sediment lost from the grid this iteration (depth in m)
Definition simulation.h:874
double m_dThisIterCliffCollapseErosionSandUncons
This-iteration total of sand unconsolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:931
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 m_dThisiterUnconsSandInput
Depth (m) of sand unconsolidated sediment added, at this iteration.
Definition simulation.h:976
double m_dAllCellsDeepWaterWavePeriod
Deep water wave period for all sea cells.
Definition simulation.h:760
double m_dThisIterActualPlatformErosionFineCons
Total actual platform erosion (fine consolidated sediment) for this iteration (depth in m)
Definition simulation.h:835
double m_dThisIterCliffCollapseErosionCoarseCons
This-iteration total of coarse consolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:943
double m_dThisIterLeftGridUnconsSand
Total sand unconsolidated sediment lost from the grid this iteration (depth in m)
Definition simulation.h:871
int nInitGridAndCalcStillWaterLevel(void)
At the beginning of each timestep: clear vector coasts, profiles, and polygons, initialize the raster...
Definition init_grid.cpp:46
double m_dStartIterSuspFineAllCells
Depth (m) of fine suspended sediment at the start of the simulation, all cells (both inside and outsi...
Definition simulation.h:982
vector< bool > m_bConsChangedThisIter
One element per layer: has the consolidated sediment of this layer been changed during this iteration...
vector< CGeomCoastPolygon * > m_pVCoastPolygon
Pointers to coast polygon objects, in down-coast sequence TODO 044 Will need to use global polygon ID...
double m_dThisIterLeftGridUnconsFine
Total fine unconsolidated sediment lost from the grid this iteration (depth in m)
Definition simulation.h:868
This file contains global definitions for CoastalME.
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:491
string const WARN
Definition cme.h:903
int const RTN_OK
Definition cme.h:694
Contains CRWCoast definitions.
Contains CGeomRasterGrid definitions.
Contains CSimulation definitions.