Skip to content

Commit 0c67b82

Browse files
MFratersMFraters
committed
Add oceanic plate euler pole velocity plugins.
1 parent fb9a7b5 commit 0c67b82

File tree

2 files changed

+241
-0
lines changed
  • include/world_builder/features/oceanic_plate_models/velocity
  • source/world_builder/features/oceanic_plate_models/velocity

2 files changed

+241
-0
lines changed

‎include/world_builder/features/oceanic_plate_models/velocity/euler_pole.h‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
/*
2+
Copyright (C) 2018-2024 by the authors of the World Builder code.
3+
4+
This file is part of the World Builder.
5+
6+
This program is free software: you can redistribute it and/or modify
7+
it under the terms of the GNU Lesser General Public License as published
8+
by the Free Software Foundation, either version 2 of the License, or
9+
(at your option) any later version.
10+
11+
This program is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU Lesser General Public License for more details.
15+
16+
You should have received a copy of the GNU Lesser General Public License
17+
along with this program. If not, see <https://www.gnu.org/licenses/>.
18+
*/
19+
20+
#ifndef WORLD_BUILDER_FEATURES_OCEANIC_PLATE_MODELS_VELOCITY_EULER_POLE_H
21+
#define WORLD_BUILDER_FEATURES_OCEANIC_PLATE_MODELS_VELOCITY_EULER_POLE_H
22+
23+
24+
#include "world_builder/features/oceanic_plate_models/velocity/interface.h"
25+
#include "world_builder/features/feature_utilities.h"
26+
#include "world_builder/objects/surface.h"
27+
28+
29+
namespace WorldBuilder
30+
{
31+
namespace Features
32+
{
33+
using namespace FeatureUtilities;
34+
namespace OceanicPlateModels
35+
{
36+
namespace Velocity
37+
{
38+
/**
39+
* This class represents a oceanic plate and can implement submodules
40+
* for velocity and composition. These submodules determine what
41+
* the returned velocity or composition of the velocity and composition
42+
* functions of this class will be.
43+
*/
44+
class EulerPole final: public Interface
45+
{
46+
public:
47+
/**
48+
* constructor
49+
*/
50+
EulerPole(WorldBuilder::World *world);
51+
52+
/**
53+
* Destructor
54+
*/
55+
~EulerPole() override final;
56+
57+
/**
58+
* declare and read in the world builder file into the parameters class
59+
*/
60+
static
61+
void declare_entries(Parameters &prm, const std::string &parent_name = "");
62+
63+
/**
64+
* declare and read in the world builder file into the parameters class
65+
*/
66+
void parse_entries(Parameters &prm, const std::vector<Point<2>> &coordinates) override final;
67+
68+
69+
/**
70+
* Returns a velocity based on the given position, depth in the model,
71+
* gravity and current velocity.
72+
*/
73+
std::array<double,3> get_velocity(const Point<3> &position,
74+
const Objects::NaturalCoordinate &position_in_natural_coordinates,
75+
const double depth,
76+
const double gravity,
77+
std::array<double,3> velocity,
78+
const double feature_min_depth,
79+
const double feature_max_depth) const override final;
80+
81+
82+
private:
83+
// euler pole velocity submodule parameters
84+
double min_depth;
85+
Objects::Surface min_depth_surface;
86+
double max_depth;
87+
Objects::Surface max_depth_surface;
88+
Point<3> euler_pole;
89+
Operations operation;
90+
91+
};
92+
} // namespace Velocity
93+
} // namespace OceanicPlateModels
94+
} // namespace Features
95+
} // namespace WorldBuilder
96+
97+
#endif

‎source/world_builder/features/oceanic_plate_models/velocity/euler_pole.cc‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
/*
2+
Copyright (C) 2018-2024 by the authors of the World Builder code.
3+
4+
This file is part of the World Builder.
5+
6+
This program is free software: you can redistribute it and/or modify
7+
it under the terms of the GNU Lesser General Public License as published
8+
by the Free Software Foundation, either version 2 of the License, or
9+
(at your option) any later version.
10+
11+
This program is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU Lesser General Public License for more details.
15+
16+
You should have received a copy of the GNU Lesser General Public License
17+
along with this program. If not, see <https://www.gnu.org/licenses/>.
18+
*/
19+
20+
#include "world_builder/features/oceanic_plate_models/velocity/euler_pole.h"
21+
22+
#include "world_builder/assert.h"
23+
#include "world_builder/consts.h"
24+
#include "world_builder/coordinate_system.h"
25+
#include "world_builder/nan.h"
26+
#include "world_builder/types/array.h"
27+
#include "world_builder/types/double.h"
28+
#include "world_builder/types/object.h"
29+
#include "world_builder/types/one_of.h"
30+
#include "world_builder/types/value_at_points.h"
31+
#include "world_builder/utilities.h"
32+
33+
34+
namespace WorldBuilder
35+
{
36+
37+
using namespace Utilities;
38+
39+
namespace Features
40+
{
41+
namespace OceanicPlateModels
42+
{
43+
namespace Velocity
44+
{
45+
EulerPole::EulerPole(WorldBuilder::World *world_)
46+
:
47+
min_depth(NaN::DSNAN),
48+
max_depth(NaN::DSNAN),
49+
euler_pole(NaN::DSNAN,NaN::DSNAN,NaN::DSNAN,CoordinateSystem::invalid),
50+
operation(Operations::REPLACE)
51+
{
52+
this->world = world_;
53+
this->name = "euler pole";
54+
}
55+
56+
EulerPole::~EulerPole()
57+
= default;
58+
59+
void
60+
EulerPole::declare_entries(Parameters &prm, const std::string & /*unused*/)
61+
{
62+
// Document plugin and require entries if needed.
63+
// Add `velocity` and to the required parameters.
64+
prm.declare_entry("", Types::Object({"euler pole"}),
65+
"Uniform velocity model. Set the velocity to a constant value.");
66+
67+
// Declare entries of this plugin
68+
prm.declare_entry("min depth", Types::OneOf(Types::Double(0),
69+
Types::Array(Types::ValueAtPoints(0.,2)),
70+
Types::String("")),
71+
"The depth in meters from which the composition of this feature is present.");
72+
73+
prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits<double>::max()),
74+
Types::Array(Types::ValueAtPoints(std::numeric_limits<double>::max(),2)),
75+
Types::String("")),
76+
"The depth in meters to which the composition of this feature is present.");
77+
78+
79+
prm.declare_entry("euler pole", Types::Array(Types::Double(0.0),2,2),
80+
"The euler pole for the plate (longitude, lattitude) in degree.");
81+
82+
prm.declare_entry("angular velocity", Types::Double(0.0), "The angular velocity of the Euler pole in degree/Myr.");
83+
84+
}
85+
86+
void
87+
EulerPole::parse_entries(Parameters &prm, const std::vector<Point<2>> &coordinates)
88+
{
89+
constexpr double year_in_seconds = 60*60*24*365.2425;
90+
WBAssertThrow(prm.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical,
91+
"The Euler pole velocity model can only be used in a spherical coordinate system.");
92+
min_depth_surface = Objects::Surface(prm.get("min depth",coordinates));
93+
min_depth = min_depth_surface.minimum;
94+
max_depth_surface = Objects::Surface(prm.get("max depth",coordinates));
95+
max_depth = max_depth_surface.maximum;
96+
operation = string_operations_to_enum(prm.get<std::string>("operation"));
97+
std::vector<double> euler_pole_degree = prm.get_vector<double>("euler pole");
98+
const double angular_speed = (1e-6/year_in_seconds)*sin(prm.get<double>("angular velocity")*Consts::PI/180.);
99+
const double longitude = euler_pole_degree[0]*Consts::PI/180.;
100+
const double lattitude = euler_pole_degree[1]*Consts::PI/180.;
101+
euler_pole[0] = angular_speed*cos(lattitude)*cos(longitude);
102+
euler_pole[1] = angular_speed*cos(lattitude)*sin(longitude);
103+
euler_pole[2] = angular_speed*sin(lattitude);
104+
}
105+
106+
107+
std::array<double,3>
108+
EulerPole::get_velocity(const Point<3> &position_in_cartesian_coordinates,
109+
const Objects::NaturalCoordinate &position_in_natural_coordinates,
110+
const double depth,
111+
const double /*gravity*/,
112+
std::array<double,3> velocity_,
113+
const double /*feature_min_depth*/,
114+
const double /*feature_max_depth*/) const
115+
{
116+
117+
if (depth <= max_depth && depth >= min_depth)
118+
{
119+
const double min_depth_local = min_depth_surface.constant_value ? min_depth : min_depth_surface.local_value(position_in_natural_coordinates.get_surface_point()).interpolated_value;
120+
const double max_depth_local = max_depth_surface.constant_value ? max_depth : max_depth_surface.local_value(position_in_natural_coordinates.get_surface_point()).interpolated_value;
121+
if (depth <= max_depth_local && depth >= min_depth_local)
122+
{
123+
const double position_norm = position_in_cartesian_coordinates.norm();
124+
const Point<3> position_normalized = position_in_cartesian_coordinates/position_norm;
125+
const Point<3> euler_velocity = position_norm * Utilities::cross_product(euler_pole, position_normalized);
126+
//std::terminate();
127+
return {{
128+
apply_operation(operation,velocity_[0],euler_velocity[0]),
129+
apply_operation(operation,velocity_[1],euler_velocity[1]),
130+
apply_operation(operation,velocity_[2],euler_velocity[2])
131+
}
132+
};
133+
}
134+
}
135+
136+
return velocity_;
137+
}
138+
139+
WB_REGISTER_FEATURE_OCEANIC_PLATE_VELOCITY_MODEL(EulerPole, euler pole)
140+
} // namespace Velocity
141+
} // namespace OceanicPlateModels
142+
} // namespace Features
143+
} // namespace WorldBuilder
144+

0 commit comments

Comments
 (0)