netxsimdg
Loading...
Searching...
No Matches
ParaGrid_test.cpp
Go to the documentation of this file.
1
8#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
9#include <doctest/doctest.h>
10
11#include "include/Configurator.hpp"
12#include "include/ConfiguredModule.hpp"
13#include "include/NZLevels.hpp"
14#include "include/ParaGridIO.hpp"
15#include "include/ParametricGrid.hpp"
16#include "include/gridNames.hpp"
17#include "include/IStructureModule.hpp"
18
19#include <cmath>
20#include <filesystem>
21#include <fstream>
22#include <sstream>
23
24#include <ncAtt.h>
25#include <ncFile.h>
26#include <ncGroup.h>
27#include <ncVar.h>
28
29const std::string filename = "paraGrid_test.nc";
30const std::string diagFile = "paraGrid_diag.nc";
31
32static const int DG = 3;
33static const int DGSTRESS = 6;
34static const int CG = 2;
35
36namespace Nextsim {
37
38size_t c = 0;
39
40TEST_SUITE_BEGIN("ParaGrid");
41TEST_CASE("Write and read a ModelState-based ParaGrid restart file")
42{
43 Module::setImplementation<IStructure>("ParametricGrid");
44
45 std::filesystem::remove(filename);
46
47 ParametricGrid grid;
48 ParaGridIO* pio = new ParaGridIO(grid);
49 grid.setIO(pio);
50
51 // Set the dimension lengths
52 size_t nx = 25;
53 size_t ny = 15;
54 size_t nz = 3;
55 NZLevels::set(nz);
56 size_t nxcg = CG * nx + 1;
57 size_t nycg = CG * ny + 1;
58
59 double yFactor = 0.01;
60 double xFactor = 0.0001;
61
62 ModelArray::setDimension(ModelArray::Dimension::X, nx);
63 ModelArray::setDimension(ModelArray::Dimension::Y, ny);
64 ModelArray::setDimension(ModelArray::Dimension::Z, NZLevels::get());
65 ModelArray::setDimension(ModelArray::Dimension::XVERTEX, nx + 1);
66 ModelArray::setDimension(ModelArray::Dimension::YVERTEX, ny + 1);
67 ModelArray::setDimension(ModelArray::Dimension::XCG, nxcg);
68 ModelArray::setDimension(ModelArray::Dimension::YCG, nycg);
69
70 ModelArray::setNComponents(ModelArray::Type::DG, DG);
71 ModelArray::setNComponents(ModelArray::Type::DGSTRESS, DGSTRESS);
72 ModelArray::setNComponents(ModelArray::Type::VERTEX, ModelArray::nCoords);
73
74 HField fractional(ModelArray::Type::H);
75 DGField fractionalDG(ModelArray::Type::DG);
76 HField mask(ModelArray::Type::H);
77 fractional.resize();
78 fractionalDG.resize();
79 for (size_t j = 0; j < ny; ++j) {
80 for (size_t i = 0; i < nx; ++i) {
81 fractional(i, j) = j * yFactor + i * xFactor;
82 mask(i, j)
83 = (i - nx / 2) * (i - nx / 2) + (j - ny / 2) * (j - ny / 2) > (nx * ny) ? 0 : 1;
84 for (size_t d = 0; d < DG; ++d) {
85 fractionalDG.components({ i, j })[d] = fractional(i, j) + d;
86 }
87 }
88 }
89
90 DGField hice = fractionalDG + 10;
91 DGField cice = fractionalDG + 20;
92 HField hsnow = fractional + 30;
93 ZField tice(ModelArray::Type::Z);
94 tice.resize();
95 for (size_t i = 0; i < ModelArray::size(ModelArray::Type::H); ++i) {
96 for (size_t k = 0; k < nz; ++k) {
97 tice.zIndexAndLayer(i, k) = fractional[i] + 40 + k;
98 }
99 }
100
101 VertexField coordinates(ModelArray::Type::VERTEX);
102 coordinates.resize();
103 // Planar coordinates
104 double scale = 1e5;
105
106 for (size_t i = 0; i < ModelArray::definedDimensions.at(ModelArray::Dimension::XVERTEX).length;
107 ++i) {
108 for (size_t j = 0;
109 j < ModelArray::definedDimensions.at(ModelArray::Dimension::YVERTEX).length; ++j) {
110 double x = i - 0.5 - nx / 2;
111 double y = j - 0.5 - ny / 2;
112 coordinates.components({ i, j })[0] = x * scale;
113 coordinates.components({ i, j })[1] = y * scale;
114 }
115 }
116
117 REQUIRE(coordinates.components({ 12, 13 })[0] - coordinates.components({ 11, 13 })[0] == scale);
118 REQUIRE(coordinates.components({ 12, 13 })[1] - coordinates.components({ 12, 12 })[1] == scale);
119
120 ModelState state = { {
121 { maskName, mask },
122 { hiceName, hice },
123 { ciceName, cice },
124 { hsnowName, hsnow },
125 { ticeName, tice },
126 { coordsName, coordinates },
127 },
128 {} };
129
130 ModelMetadata metadata;
131 metadata.setTime(TimePoint("2000-01-01T00:00:00Z"));
132
133 grid.dumpModelState(state, metadata, filename, true);
134
135 REQUIRE(std::filesystem::exists(std::filesystem::path(filename)));
136
137 // Reset the array dimensions to make sure that the read function gets them correct
138 ModelArray::setDimension(ModelArray::Dimension::X, 1);
139 ModelArray::setDimension(ModelArray::Dimension::Y, 1);
140 ModelArray::setDimension(ModelArray::Dimension::Z, 1);
141 ModelArray::setDimension(ModelArray::Dimension::XVERTEX, 1);
142 ModelArray::setDimension(ModelArray::Dimension::YVERTEX, 1);
143 ModelArray::setDimension(ModelArray::Dimension::XCG, 1);
144 ModelArray::setDimension(ModelArray::Dimension::YCG, 1);
145 // In the full model numbers of DG components are set at compile time, so they are not reset
146 REQUIRE(ModelArray::nComponents(ModelArray::Type::DG) == DG);
147 REQUIRE(ModelArray::nComponents(ModelArray::Type::VERTEX) == ModelArray::nCoords);
148
149 ParametricGrid gridIn;
150 ParaGridIO* readIO = new ParaGridIO(gridIn);
151 gridIn.setIO(readIO);
152
153 ModelState ms = gridIn.getModelState(filename);
154
155 REQUIRE(ModelArray::dimensions(ModelArray::Type::Z)[0] == nx);
156 REQUIRE(ModelArray::dimensions(ModelArray::Type::Z)[1] == ny);
157 REQUIRE(ModelArray::dimensions(ModelArray::Type::Z)[2] == NZLevels::get());
158
159 REQUIRE(ms.data.size() == state.data.size());
160
161 ModelArray& ticeRef = ms.data.at(ticeName);
162 REQUIRE(ModelArray::nDimensions(ModelArray::Type::Z) == 3);
163 REQUIRE(ticeRef.getType() == ModelArray::Type::Z);
164 REQUIRE(ticeRef.nDimensions() == 3);
165 REQUIRE(ticeRef.dimensions()[0] == nx);
166 REQUIRE(ticeRef.dimensions()[1] == ny);
167 REQUIRE(ticeRef.dimensions()[2] == NZLevels::get());
168
169 ModelArray& hiceRef = ms.data.at(hiceName);
170 REQUIRE(hiceRef.nDimensions() == 2);
171 REQUIRE(hiceRef.dimensions()[0] == nx);
172 REQUIRE(hiceRef.dimensions()[1] == ny);
173 REQUIRE(ModelArray::nComponents(ModelArray::Type::DG) == DG);
174 REQUIRE(hiceRef.nComponents() == DG);
175
176 REQUIRE(ticeRef(12, 14, 1) == tice(12, 14, 1));
177
178 ModelArray& coordRef = ms.data.at(coordsName);
179 REQUIRE(coordRef.nDimensions() == 2);
180 REQUIRE(coordRef.nComponents() == 2);
181 REQUIRE(coordRef.components({ 12, 13 })[0] - coordRef.components({ 11, 13 })[0] == scale);
182 REQUIRE(coordRef.components({ 12, 13 })[1] - coordRef.components({ 12, 12 })[1] == scale);
183
184 std::filesystem::remove(filename);
185}
186
187TEST_CASE("Write a diagnostic ParaGrid file")
188{
189 Module::setImplementation<IStructure>("ParametricGrid");
190
191 REQUIRE(Module::getImplementation<IStructure>().structureType() == "parametric_rectangular");
192
193 std::filesystem::remove(diagFile);
194
195 ParametricGrid grid;
196 ParaGridIO* pio = new ParaGridIO(grid);
197 grid.setIO(pio);
198
199 // Set the dimension lengths
200 size_t nx = 30;
201 size_t ny = 20;
202 size_t nz = 3;
203 NZLevels::set(nz);
204 size_t nxcg = CG * nx + 1;
205 size_t nycg = CG * ny + 1;
206
207 double yFactor = 0.01;
208 double xFactor = 0.0001;
209
210 ModelArray::setDimension(ModelArray::Dimension::X, nx);
211 ModelArray::setDimension(ModelArray::Dimension::Y, ny);
212 ModelArray::setDimension(ModelArray::Dimension::Z, NZLevels::get());
213 ModelArray::setDimension(ModelArray::Dimension::XVERTEX, nx + 1);
214 ModelArray::setDimension(ModelArray::Dimension::YVERTEX, ny + 1);
215 ModelArray::setDimension(ModelArray::Dimension::XCG, nxcg);
216 ModelArray::setDimension(ModelArray::Dimension::YCG, nycg);
217
218 ModelArray::setNComponents(ModelArray::Type::DG, DG);
219 ModelArray::setNComponents(ModelArray::Type::DGSTRESS, DGSTRESS);
220 ModelArray::setNComponents(ModelArray::Type::VERTEX, ModelArray::nCoords);
221
222 HField fractional(ModelArray::Type::H);
223 DGField fractionalDG(ModelArray::Type::DG);
224 HField mask(ModelArray::Type::H);
225 fractional.resize();
226 fractionalDG.resize();
227 for (size_t j = 0; j < ny; ++j) {
228 for (size_t i = 0; i < nx; ++i) {
229 fractional(i, j) = j * yFactor + i * xFactor;
230 mask(i, j) = fractional(i, j);
231// = (i - nx / 2) * (i - nx / 2) + (j - ny / 2) * (j - ny / 2) > (nx * ny) ? 0 : 1;
232 for (size_t d = 0; d < DG; ++d) {
233 fractionalDG.components({ i, j })[d] = fractional(i, j) + d;
234 }
235 }
236 }
237 double prec = 1e-9;
238 REQUIRE(fractional(12, 12) - fractional(11, 12) == doctest::Approx(xFactor).epsilon(prec));
239 REQUIRE(fractional(12, 12) - fractional(12, 11) == doctest::Approx(yFactor).epsilon(prec));
240
241 REQUIRE(fractionalDG(12, 12) - fractionalDG(11, 12) == doctest::Approx(xFactor).epsilon(prec));
242 REQUIRE(fractionalDG(12, 12) - fractionalDG(12, 11) == doctest::Approx(yFactor).epsilon(prec));
243
244
245 DGField hice = fractionalDG + 10;
246 DGField cice = fractionalDG + 20;
247
248 VertexField coordinates(ModelArray::Type::VERTEX);
249 coordinates.resize();
250 // Planar coordinates
251 double scale = 1e5;
252
253 for (size_t i = 0; i < ModelArray::definedDimensions.at(ModelArray::Dimension::XVERTEX).length;
254 ++i) {
255 for (size_t j = 0;
256 j < ModelArray::definedDimensions.at(ModelArray::Dimension::YVERTEX).length; ++j) {
257 double x = i - 0.5 - nx / 2;
258 double y = j - 0.5 - ny / 2;
259 coordinates.components({ i, j })[0] = x * scale;
260 coordinates.components({ i, j })[1] = y * scale;
261 }
262 }
263
264 REQUIRE(coordinates.components({ 12, 13 })[0] - coordinates.components({ 11, 13 })[0] == scale);
265 REQUIRE(coordinates.components({ 12, 13 })[1] - coordinates.components({ 12, 12 })[1] == scale);
266
267 ModelState state = { {
268 { maskName, mask },
269 { hiceName, hice },
270 { ciceName, cice },
271 { coordsName, coordinates },
272 },
273 {} };
274
275 ModelMetadata metadata;
276 metadata.setTime(TimePoint("2000-01-01T00:00:00Z"));
277
278 grid.dumpModelState(state, metadata, diagFile, false);
279
280 for (int t = 1; t < 5; ++t) {
281 hice += 100;
282 cice += 100;
283 state = { {
284 { hiceName, hice },
285 { ciceName, cice },
286 },
287 {} };
288
289 grid.dumpModelState(state, metadata, diagFile, false);
290 }
291 pio->close(diagFile);
292
293 // What do we have in the file?
294 netCDF::NcFile ncFile(diagFile, netCDF::NcFile::read);
295
296 REQUIRE(ncFile.getGroups().size() == 3);
297 netCDF::NcGroup structGrp(ncFile.getGroup(IStructure::structureNodeName()));
298 netCDF::NcGroup metaGrp(ncFile.getGroup(IStructure::metadataNodeName()));
299 netCDF::NcGroup dataGrp(ncFile.getGroup(IStructure::dataNodeName()));
300
301 std::string structureType;
302 structGrp.getAtt(grid.typeNodeName()).getValues(structureType);
303 // TODO implement ParametricGrid in the IStructureModule
304 REQUIRE(structureType == grid.structureType());
305
306 // TODO test metadata
307
308 // test data
309 REQUIRE(dataGrp.getVarCount() == 5);
310 netCDF::NcVar hiceVar = dataGrp.getVar(hiceName);
311 netCDF::NcVar ciceVar = dataGrp.getVar(ciceName);
312 netCDF::NcVar coordVar = dataGrp.getVar(coordsName);
313 netCDF::NcVar maskVar = dataGrp.getVar(maskName);
314 netCDF::NcVar timeVar = dataGrp.getVar(timeName);
315
316 // hice
317 REQUIRE(hiceVar.getDimCount() == 4);
318
319 ncFile.close();
320
321 std::filesystem::remove(diagFile);
322
323}
324
325#define TO_STR(s) TO_STRI(s)
326#define TO_STRI(s) #s
327#ifndef TEST_FILE_SOURCE
328#define TEST_FILE_SOURCE .
329#endif
330
331TEST_CASE("Test array ordering")
332{
333 std::string inputFilename = "ParaGridIO_input_test.nc";
334
335 Module::setImplementation<IStructure>("ParametricGrid");
336
337 REQUIRE(Module::getImplementation<IStructure>().structureType() == "parametric_rectangular");
338
339 size_t nx = 9;
340 size_t ny = 11;
341 NZLevels::set(1);
342
343 double xFactor = 10;
344
345 ModelArray::setDimension(ModelArray::Dimension::X, nx);
346 ModelArray::setDimension(ModelArray::Dimension::Y, ny);
347 ModelArray::setDimension(ModelArray::Dimension::Z, NZLevels::get());
348
349 HField index2d(ModelArray::Type::H);
350 index2d.resize();
351 std::string fieldName = "index2d";
352 std::set<std::string> fields = { fieldName };
353 TimePoint time;
354
355 ModelState state = ParaGridIO::readForcingTimeStatic(fields, time, TO_STR(TEST_FILE_SOURCE) + std::string("/") + inputFilename);
356 REQUIRE(state.data.count(fieldName) > 0);
357 index2d = state.data.at(fieldName);
358 REQUIRE(index2d(3, 5) == 35);
359 // And that's all that's needed
360}
361
362#undef TO_STR
363#undef TO_STRI
364TEST_SUITE_END();
365
366}
Nextsim::IStructure::dataNodeName
static const std::string dataNodeName()
Returns the name of the data node.
Definition IStructure.hpp:84
Nextsim::IStructure::structureNodeName
static const std::string structureNodeName()
The name of the group holding the definitive structure type.
Definition IStructure.hpp:86
Nextsim::IStructure::metadataNodeName
static const std::string metadataNodeName()
Returns the name of the metadata node.
Definition IStructure.hpp:82
Nextsim::ModelArray::dimensions
const MultiDim & dimensions() const
Returns a vector<size_t> of the size of each dimension of this type of ModelArray.
Definition ModelArray.hpp:225
Nextsim::ModelArray::setDimension
static void setDimension(Dimension dim, size_t length)
Sets the length of an individual dimension before propagating it to the defined array types.
Definition ModelArray.cpp:203
Nextsim::ModelArray::nComponents
size_t nComponents() const
Returns the number of discontinuous Galerkin components held in this type of ModelArray.
Definition ModelArray.hpp:532
Nextsim::ModelArray::nDimensions
size_t nDimensions() const
Returns the number of dimensions of this type of ModelArray.
Definition ModelArray.hpp:221
Nextsim::ModelArray::size
size_t size() const
Returns the total number of elements of this type of ModelArray.
Definition ModelArray.hpp:229
Nextsim::ModelArray::definedDimensions
static std::map< Dimension, DimensionSpec > definedDimensions
The name and length of each dimension that is defined.
Definition ModelArray.hpp:66
Nextsim::ModelArray::setNComponents
static void setNComponents(std::map< Type, size_t > cMap)
Sets the number of components for DG & CG array types.
Definition ModelArray.cpp:196
Ice::epsilon
const double epsilon
Thermal emissivity of smooth ice [0..1].
Definition constants.hpp:43