Skip to content

Commit 91e0867

Browse files
MFratersMFraters
committed
WIP
1 parent 5cc6fa2 commit 91e0867

File tree

9 files changed

+111
-44
lines changed

9 files changed

+111
-44
lines changed

‎source/gwb-grid/main.cc‎

Lines changed: 52 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
#include "world_builder/assert.h"
2929
#include "world_builder/coordinate_system.h"
3030
#include "world_builder/nan.h"
31+
#include "world_builder/objects/natural_coordinate.h"
3132
#include "world_builder/point.h"
3233
#include "world_builder/utilities.h"
3334
#include "world_builder/world.h"
@@ -921,13 +922,13 @@ int main(int argc, char **argv)
921922

922923
const double dlong = opening_angle_long_rad / static_cast<double>(n_cell_x);
923924
const double dlat = dim == 3 ? opening_angle_lat_rad / static_cast<double>(n_cell_y) : 0.;
924-
const double lr = outer_radius - inner_radius;
925-
const double dr = lr / static_cast<double>(n_cell_z);
926925

927926
grid_x.resize(n_p);
928927
grid_y.resize(dim == 3 ? n_p : 0);
929928
grid_z.resize(n_p);
930929
grid_depth.resize(n_p);
930+
std::vector<double> domain_height(n_p);
931+
std::vector<double> cell_height(n_p);
931932

932933
std::cout << "[4/6] Building the grid: stage 1 of 3 \r";
933934
std::cout.flush();
@@ -938,8 +939,10 @@ int main(int argc, char **argv)
938939
for (size_t j = 1; j <= n_cell_z + 1; ++j)
939940
{
940941
grid_x[counter] = x_min + (static_cast<double>(i) - 1.0) * dlong;
941-
grid_z[counter] = inner_radius + (static_cast<double>(j) - 1.0) * dr;
942-
grid_depth[counter] = lr - (static_cast<double>(j) - 1.0) * dr;
942+
domain_height [counter]= outer_radius - inner_radius;
943+
cell_height[counter] = domain_height[counter] / static_cast<double>(n_cell_z);
944+
grid_z[counter] = inner_radius + (static_cast<double>(j) - 1.0) * cell_height[counter];
945+
grid_depth[counter] = domain_height[counter] - (static_cast<double>(j) - 1.0) * cell_height[counter];
943946
counter++;
944947
}
945948
}
@@ -953,8 +956,16 @@ int main(int argc, char **argv)
953956
{
954957
grid_x[counter] = x_min + (static_cast<double>(i) - 1.0) * dlong;
955958
grid_y[counter] = y_min + (static_cast<double>(j) - 1.0) * dlat;
956-
grid_z[counter] = inner_radius + (static_cast<double>(k) - 1.0) * dr;
957-
grid_depth[counter] = lr - (static_cast<double>(k) - 1.0) * dr;
959+
// todo: actual position if that is what we decide on?
960+
std::vector<std::array<unsigned ,3>> properties;
961+
properties.push_back({{6,0,0}}); // topography
962+
const std::array<double,3> coords = {{grid_x[counter], grid_y[counter],inner_radius}};
963+
world->properties(coords, 0,properties);
964+
const double topography = world->properties(coords, 0,properties)[0];
965+
domain_height[counter] = outer_radius + topography - inner_radius;
966+
cell_height[counter] = domain_height[counter] / static_cast<double>(n_cell_z);
967+
grid_z[counter] = inner_radius + (static_cast<double>(k) - 1.0) * cell_height[counter];
968+
grid_depth[counter] = domain_height[counter] - (static_cast<double>(k) - 1.0) * cell_height[counter];
958969
counter++;
959970
}
960971
}
@@ -970,50 +981,50 @@ int main(int argc, char **argv)
970981
// position 0 of this cell
971982
grid_x[counter] = x_min + static_cast<double>(i) * dlong;
972983
grid_y[counter] = y_min + static_cast<double>(j) * dlat;
973-
grid_z[counter] = inner_radius + static_cast<double>(k) * dr;
974-
grid_depth[counter] = lr - static_cast<double>(k) * dr;
984+
grid_z[counter] = inner_radius + static_cast<double>(k) * cell_height[counter];
985+
grid_depth[counter] = domain_height[counter] - static_cast<double>(k) * cell_height[counter];
975986
counter++;
976987
// position 1 of this cell
977988
grid_x[counter] = x_min + (static_cast<double>(i) + 1.0) * dlong;
978989
grid_y[counter] = y_min + static_cast<double>(j) * dlat;
979-
grid_z[counter] = inner_radius + static_cast<double>(k) * dr;
980-
grid_depth[counter] = lr - static_cast<double>(k) * dr;
990+
grid_z[counter] = inner_radius + static_cast<double>(k) * cell_height[counter];
991+
grid_depth[counter] = domain_height[counter] - static_cast<double>(k) * cell_height[counter];
981992
counter++;
982993
// position 2 of this cell
983994
grid_x[counter] = x_min + (static_cast<double>(i) + 1.0) * dlong;
984995
grid_y[counter] = y_min + (static_cast<double>(j) + 1.0) * dlat;
985-
grid_z[counter] = inner_radius + static_cast<double>(k) * dr;
986-
grid_depth[counter] = lr - static_cast<double>(k) * dr;
996+
grid_z[counter] = inner_radius + static_cast<double>(k) * cell_height[counter];
997+
grid_depth[counter] = domain_height[counter] - static_cast<double>(k) * cell_height[counter];
987998
counter++;
988999
// position 3 of this cell
9891000
grid_x[counter] = x_min + static_cast<double>(i) * dlong;
9901001
grid_y[counter] = y_min + (static_cast<double>(j) + 1.0) * dlat;
991-
grid_z[counter] = inner_radius + static_cast<double>(k) * dr;
992-
grid_depth[counter] = lr - static_cast<double>(k) * dr;
1002+
grid_z[counter] = inner_radius + static_cast<double>(k) * cell_height[counter];
1003+
grid_depth[counter] = domain_height[counter] - static_cast<double>(k) * cell_height[counter];
9931004
counter++;
9941005
// position 0 of this cell
9951006
grid_x[counter] = x_min + static_cast<double>(i) * dlong;
9961007
grid_y[counter] = y_min + static_cast<double>(j) * dlat;
997-
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * dr;
998-
grid_depth[counter] = lr - (static_cast<double>(k) + 1.0) * dr;
1008+
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * cell_height[counter];
1009+
grid_depth[counter] = domain_height[counter] - (static_cast<double>(k) + 1.0) * cell_height[counter];
9991010
counter++;
10001011
// position 1 of this cell
10011012
grid_x[counter] = x_min + (static_cast<double>(i) + 1.0) * dlong;
10021013
grid_y[counter] = y_min + static_cast<double>(j) * dlat;
1003-
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * dr;
1004-
grid_depth[counter] = lr - (static_cast<double>(k) + 1.0) * dr;
1014+
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * cell_height[counter];
1015+
grid_depth[counter] = domain_height[counter] - (static_cast<double>(k) + 1.0) * cell_height[counter];
10051016
counter++;
10061017
// position 2 of this cell
10071018
grid_x[counter] = x_min + (static_cast<double>(i) + 1.0) * dlong;
10081019
grid_y[counter] = y_min + (static_cast<double>(j) + 1.0) * dlat;
1009-
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * dr;
1010-
grid_depth[counter] = lr - (static_cast<double>(k) + 1.0) * dr;
1020+
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * cell_height[counter];
1021+
grid_depth[counter] = domain_height[counter] - (static_cast<double>(k) + 1.0) * cell_height[counter];
10111022
counter++;
10121023
// position 3 of this cell
10131024
grid_x[counter] = x_min + static_cast<double>(i) * dlong;
10141025
grid_y[counter] = y_min + (static_cast<double>(j) + 1.0) * dlat;
1015-
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * dr;
1016-
grid_depth[counter] = lr - (static_cast<double>(k) + 1.0) * dr;
1026+
grid_z[counter] = inner_radius + (static_cast<double>(k) + 1.0) * cell_height[counter];
1027+
grid_depth[counter] = domain_height[counter] - (static_cast<double>(k) + 1.0) * cell_height[counter];
10171028
WBAssert(counter < n_p, "Assert counter smaller then n_P: counter = " << counter << ", n_p = " << n_p);
10181029
counter++;
10191030
}
@@ -1567,6 +1578,7 @@ int main(int argc, char **argv)
15671578
std::vector<vtu11::DataSetInfo> dataSetInfo
15681579
{
15691580
{ "Depth", vtu11::DataSetType::PointData, 1 },
1581+
{ "Topography", vtu11::DataSetType::PointData, 1 },
15701582
{ "Temperature", vtu11::DataSetType::PointData, 1 },
15711583
{ "velocity", vtu11::DataSetType::PointData, 3 },
15721584
{ "Tag", vtu11::DataSetType::PointData, 1 },
@@ -1580,6 +1592,7 @@ int main(int argc, char **argv)
15801592
std::cout.flush();
15811593

15821594
std::vector<std::array<unsigned ,3>> properties;
1595+
properties.push_back({{6,0,0}}); // topography
15831596
properties.push_back({{1,0,0}}); // temperature
15841597

15851598
properties.push_back({{5,0,0}}); // velocity
@@ -1591,13 +1604,14 @@ int main(int argc, char **argv)
15911604

15921605

15931606
// compute temperature
1594-
std::vector<vtu11::DataSetData> data_set(4+compositions);
1607+
std::vector<vtu11::DataSetData> data_set(5+compositions);
15951608
data_set[0] = grid_depth;
15961609
data_set[1].resize(n_p);
1597-
data_set[2].resize(n_p*3);
1598-
data_set[3].resize(n_p);
1610+
data_set[2].resize(n_p);
1611+
data_set[3].resize(n_p*3);
1612+
data_set[4].resize(n_p);
15991613
for (size_t c = 0; c < compositions; ++c)
1600-
data_set[4+c].resize(n_p);
1614+
data_set[5+c].resize(n_p);
16011615

16021616
if (dim == 2)
16031617
{
@@ -1606,13 +1620,14 @@ int main(int argc, char **argv)
16061620
const std::array<double,2> coords = {{grid_x[i], grid_z[i]}};
16071621
std::vector<double> output = world->properties(coords, grid_depth[i],properties);
16081622
data_set[1][i] = output[0];
1609-
data_set[2][3*i] = output[1];
1610-
data_set[2][3*i+1] = output[2];
1611-
data_set[2][3*i+2] = output[3];
1612-
data_set[3][i] = output[4];
1623+
data_set[2][i] = output[1];
1624+
data_set[3][3*i] = output[2];
1625+
data_set[3][3*i+1] = output[3];
1626+
data_set[3][3*i+2] = output[4];
1627+
data_set[4][i] = output[5];
16131628
for (size_t c = 0; c < compositions; ++c)
16141629
{
1615-
data_set[4+c][i] = output[5+c];
1630+
data_set[5+c][i] = output[6+c];
16161631
}
16171632
});
16181633
}
@@ -1623,13 +1638,14 @@ int main(int argc, char **argv)
16231638
const std::array<double,3> coords = {{grid_x[i], grid_y[i], grid_z[i]}};
16241639
std::vector<double> output = world->properties(coords, grid_depth[i],properties);
16251640
data_set[1][i] = output[0];
1626-
data_set[2][3*i] = output[1];
1627-
data_set[2][3*i+1] = output[2];
1628-
data_set[2][3*i+2] = output[3];
1629-
data_set[3][i] = output[4];
1641+
data_set[2][i] = output[1];
1642+
data_set[3][3*i] = output[2];
1643+
data_set[3][3*i+1] = output[3];
1644+
data_set[3][3*i+2] = output[4];
1645+
data_set[4][i] = output[5];
16301646
for (size_t c = 0; c < compositions; ++c)
16311647
{
1632-
data_set[4+c][i] = output[5+c];
1648+
data_set[5+c][i] = output[6+c];
16331649
}
16341650
});
16351651
}

‎source/world_builder/features/continental_plate.cc‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -295,6 +295,11 @@ namespace WorldBuilder
295295
//std::cout << "vel=" << output[entry_in_output[i_property]] << ":" << output[entry_in_output[i_property]+1] << ":" << output[entry_in_output[i_property]+2] << std::endl;
296296
break;
297297
}
298+
case 6: // topography
299+
{
300+
// the topography is not handled by the individual features.
301+
break;
302+
}
298303
default:
299304
{
300305
WBAssertThrow(false,

‎source/world_builder/features/fault.cc‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -763,6 +763,11 @@ namespace WorldBuilder
763763
output[entry_in_output[i_property]+2] = velocity_current_section[2] + section_fraction * (velocity_next_section[2] - velocity_current_section[2]);
764764
break;
765765
}
766+
case 6: // topography
767+
{
768+
// the topography is not handled by the individual features.
769+
break;
770+
}
766771
default:
767772
{
768773
WBAssertThrow(false,

‎source/world_builder/features/mantle_layer.cc‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -288,7 +288,12 @@ namespace WorldBuilder
288288
output[entry_in_output[i_property]+2] = velocity[2];
289289
break;
290290
}
291-
break;
291+
case 6: // topography
292+
{
293+
// the topography is not handled by the individual features.
294+
break;
295+
}
296+
292297
default:
293298
{
294299
WBAssertThrow(false,

‎source/world_builder/features/oceanic_plate.cc‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -306,7 +306,12 @@ namespace WorldBuilder
306306
output[entry_in_output[i_property]+2] = velocity[2];
307307
break;
308308
}
309-
break;
309+
case 6: // topography
310+
{
311+
// the topography is not handled by the individual features.
312+
break;
313+
}
314+
310315
default:
311316
{
312317
WBAssertThrow(false,

‎source/world_builder/features/plume.cc‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -428,6 +428,11 @@ namespace WorldBuilder
428428
output[entry_in_output[i_property]+2] = velocity[2];
429429
break;
430430
}
431+
case 6: // topography
432+
{
433+
// the topography is not handled by the individual features.
434+
break;
435+
}
431436
default:
432437
{
433438
WBAssertThrow(false,
@@ -443,4 +448,4 @@ namespace WorldBuilder
443448
WB_REGISTER_FEATURE(Plume, plume)
444449

445450
} // namespace Features
446-
} // namespace WorldBuilder
451+
} // namespace WorldBuilder

‎source/world_builder/features/subducting_plate.cc‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -796,6 +796,11 @@ namespace WorldBuilder
796796
output[entry_in_output[i_property]+2] = velocity_current_section[2] + section_fraction * (velocity_next_section[2] - velocity_current_section[2]);
797797
break;
798798
}
799+
case 6: // topography
800+
{
801+
// the topography is not handled by the individual features.
802+
break;
803+
}
799804
default:
800805
{
801806
WBAssertThrow(false,

‎source/world_builder/topography/uniform.cc‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,8 +38,12 @@ namespace WorldBuilder
3838
= default;
3939

4040
void
41-
Uniform::declare_entries(Parameters &/*prm*/, const std::string & /*unused*/)
41+
Uniform::declare_entries(Parameters &prm, const std::string & /*unused*/)
4242
{
43+
prm.declare_entry("", Types::Object(),
44+
"Uniform gravity model. It returns the gravity vector in a Cartesian coordinate system at "
45+
"a given position, which has a constant magitude for the whole domain. The vector points "
46+
"down in cartesian coordinates and to the center of the sphere in spherical coordinates.");
4347
// Nothing to declare.
4448
}
4549

@@ -58,7 +62,7 @@ namespace WorldBuilder
5862
const Objects::NaturalCoordinate &/*position_in_natural_coordinates*/,
5963
const double /*depth*/) const
6064
{
61-
return 0;
65+
return 10e3;
6266
}
6367

6468
Uniform::SurfaceType

‎source/world_builder/world.cc‎

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -313,10 +313,15 @@ namespace WorldBuilder
313313
n_output_entries += 3;
314314
break;
315315
}
316+
case 6: // topography (1 entry)
317+
{
318+
n_output_entries += 1;
319+
break;
320+
}
316321
default:
317322
WBAssertThrow(false,
318323
"Internal error: Unimplemented property provided. " <<
319-
"Only temperature (1), composition (2), grains (3), tag (4) or velocity (5) are allowed. "
324+
"Only temperature (1), composition (2), grains (3), tag (4), velocity (5) or topography are allowed. "
320325
"Provided property number was: " << property[0]);
321326
}
322327
}
@@ -399,11 +404,16 @@ namespace WorldBuilder
399404
counter += 3;
400405
break;
401406
}
407+
case 6: // topography
408+
{
409+
counter += 1;
410+
break;
411+
}
402412
default:
403413
{
404414
WBAssert(false,
405415
"Internal error: Unimplemented property provided by internal process. " <<
406-
"Only temperature (1), composition (2), grains (3), tag (4) or velocity (5) are allowed. "
416+
"Only temperature (1), composition (2), grains (3), tag (4), velocity (5) or topography are allowed. "
407417
"Provided property number was: " << property[0]);
408418
}
409419
}
@@ -485,10 +495,17 @@ namespace WorldBuilder
485495
properties_local.emplace_back(properties[i_property]);
486496
break;
487497
}
498+
case 6: // topography
499+
{
500+
entry_in_output.emplace_back(output.size());
501+
output.emplace_back(this->parameters.topography_model->get_topography(point, natural_coordinate, depth));
502+
properties_local.emplace_back(properties[i_property]);
503+
break;
504+
}
488505
default:
489506
WBAssertThrow(false,
490507
"Internal error: Unimplemented property provided. " <<
491-
"Only temperature (1), composition (2), grains (3), tag (4) or velocity (5) are allowed. "
508+
"Only temperature (1), composition (2), grains (3), tag (4), velocity (5) or topography (6) are allowed. "
492509
"Provided property number was: " << properties[i_property][0]);
493510
}
494511
}

0 commit comments

Comments
 (0)