Commit bcac75b6 authored by Marcus Mohr's avatar Marcus Mohr
Browse files

WIP: First draft of plate velocity computations

This draft does not, yet, compile and needs some further pimping to be
HyTeG-style. It is based on:

https://gitlab.lrz.de/marcus.mohr/plates-refactoring/-/tags/HYTEG_INCLUSION_CANDIDATE_1
parent 7e5a1344
Pipeline #40134 canceled with stages
in 23 seconds
...@@ -8,4 +8,7 @@ add_library( terraneo ) ...@@ -8,4 +8,7 @@ add_library( terraneo )
target_link_libraries( terraneo PUBLIC core hyteg ) target_link_libraries( terraneo PUBLIC core hyteg )
# target_sources( terraneo ) # target_sources( terraneo )
add_subdirectory( dataimport )
add_subdirectory( helpers )
add_subdirectory( plates )
add_subdirectory( sphericalharmonics ) add_subdirectory( sphericalharmonics )
target_sources( hyteg
PRIVATE
io.hpp
)
/*
* Copyright (c) 2022 Berta Vilacis, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program 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.
*
* 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.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include <fstream>
#include "core/extern/json.hpp"
#include "terraneo/plates/types.hpp"
namespace terraneo {
namespace io {
using walberla::real_c;
using walberla::real_t;
using walberla::uint_t;
/// Auxilliary function for opening an input file
/// \param fileName name (including potentially a path) of the input file
/// \return ifstream for the file
std::ifstream openFileForReading( std::string& fileName )
{
std::ifstream infile( fileName, std::ios::in );
if ( !infile.is_open() )
{
std::cerr << "Failed to open file '" << fileName << "' for reading!" << std::endl;
std::abort();
}
return infile;
}
/// Imports a file in JSON format
nlohmann::json readJsonFile( std::string filename )
{
std::ifstream infile{ openFileForReading( filename ) };
nlohmann::json inobj;
infile >> inobj;
return inobj;
}
/// Imports data from a text file with rotation information
std::vector< terraneo::plates::RotationInfo > readRotationsFile( std::string fileName )
{
std::ifstream file{ openFileForReading( fileName ) };
// Determine number of lines in file for reserving space
std::string line;
uint_t numRotations{ 0 };
while ( std::getline( file, line ) )
{
numRotations++;
}
file.clear(); // clear EOF
file.seekg( 0 ); // rewind stream
std::cout << "Detected " << numRotations << " rotations" << std::endl;
std::vector< terraneo::plates::RotationInfo > rotations( numRotations, terraneo::plates::RotationInfo() );
std::size_t charsUsed;
uint_t k{ 0 };
// this was primarily introduced for the regression testing during the
// refactoring of the original coding attempt; the accuracy of the data
// in the input file is smaller than even binary32 ;-)
#ifdef WALBERLA_DOUBLE_ACCURACY
#define PLATES_IO_STR_TO_FP std::stod
#else
#define PLATES_IO_STR_TO_FP std::stof
#endif
while ( std::getline( file, line ) )
{
// position inside line-string
std::size_t pos{ 0 };
rotations[k].plateID = std::stoul( line, &charsUsed );
pos += charsUsed;
rotations[k].time = real_c( PLATES_IO_STR_TO_FP( line.substr( pos ), &charsUsed ) );
pos += charsUsed;
rotations[k].latitude = real_c( PLATES_IO_STR_TO_FP( line.substr( pos ), &charsUsed ) );
pos += charsUsed;
rotations[k].longitude = real_c( PLATES_IO_STR_TO_FP( line.substr( pos ), &charsUsed ) );
pos += charsUsed;
rotations[k].angle = real_c( PLATES_IO_STR_TO_FP( line.substr( pos ), &charsUsed ) );
pos += charsUsed;
rotations[k].conjugateID = std::stoul( line.substr( pos ), &charsUsed );
// don't forget to increment index into vector
++k;
}
// Safety check
if ( k != numRotations )
{
std::cout << " Imported " << k << " rotations,\n"
<< " but should have been " << numRotations << std::endl;
std::abort();
}
return rotations;
}
} // namespace io
} // namespace terraneo
target_sources( hyteg
PRIVATE
conversions.hpp
typeAliases.hpp
)
/*
* Copyright (c) 2022 Berta Vilacis, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program 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.
*
* 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.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "core/math/Constants.h"
#include "terraneo/helpers/typeAliases.hpp"
namespace terraneo::conversions {
using walberla::math::pi;
/// Transforms angle from degrees to radians
inline real_t degToRad( real_t degree )
{
return ( degree * ( pi / real_c( 180 ) ) );
}
/// Transforms vector of angles componentwise from degrees to radians
inline vec3D degToRad( vec3D degree )
{
return ( degree * ( pi / real_c( 180 ) ) );
}
/// Transforms angle from radians to degrees
inline real_t radToDeg( real_t radian )
{
return ( radian * ( real_c( 180 ) / pi ) );
}
/// Transforms vector of angles componentwise from radians to degrees
inline vec3D radToDeg( vec3D radian )
{
return ( radian * ( real_c( 180 ) / pi ) );
}
/// Transform 3D vector from spherical to cartesian coordintates
///
/// Transform 3D vector from spherical coordinates (lon, lat, rad) to cartesian
/// ones (x,y,z), is radius is not given assume unit point on sphere.
inline vec3D sph2cart( const std::vector< real_t >& lonlat, const real_t radius = real_c( 1 ) )
{
vec3D xyz;
xyz[0] = radius * cos( degToRad( lonlat[1] ) ) * cos( degToRad( lonlat[0] ) );
xyz[1] = radius * cos( degToRad( lonlat[1] ) ) * sin( degToRad( lonlat[0] ) );
xyz[2] = radius * sin( degToRad( lonlat[1] ) );
return xyz;
}
/// Transform 3D vector from cartesian (x, y, z) to spherical coordinates (lon, lat, rad)
vec3D cart2sph( const vec3D& xyz )
{
vec3D lonlatrad;
lonlatrad[0] = radToDeg( atan2( xyz[1], xyz[0] ) );
lonlatrad[1] = radToDeg( atan2( xyz[2], sqrt( xyz[0] * xyz[0] + xyz[1] * xyz[1] ) ) );
lonlatrad[2] = xyz.norm();
return lonlatrad;
}
} // namespace terraneo::conversions
/*
* Copyright (c) 2022 Berta Vilacis, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program 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.
*
* 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.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "hyteg/eigen/EigenWrapper.hpp"
namespace terraneo {
using walberla::real_c;
using walberla::real_t;
using walberla::uint_t;
using vec3D = Eigen::Vector3d;
using mat3D = Eigen::Matrix3d;
} // namespace terraneo
target_sources( hyteg
PRIVATE
functionsForGeometry.hpp
functionsForPlates.hpp
functionsForRotations.hpp
PlateRotationProvider.hpp
PlateStorage.hpp
PlateVelocityProvider.hpp
SmoothingStrategies.hpp
types.hpp
)
/*
* Copyright (c) 2022 Berta Vilacis, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program 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.
*
* 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.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "terraneo/dataimport/io.hpp"
namespace terraneo {
namespace plates {
/// Class storing and providing access to plate rotation information
///
/// This class deals with importating infromation reconstructed plate
/// rotation information from an external datafile, stores that information
/// and provides access to it.
class PlateRotationProvider
{
public:
template < typename ImportStrategy >
PlateRotationProvider( std::string nameOfRotationsFile, ImportStrategy readRotationsFile )
{
rotInfos_ = readRotationsFile( nameOfRotationsFile );
}
const std::vector< RotationInfo >& getRotations() const { return rotInfos_; };
private:
std::vector< RotationInfo > rotInfos_;
};
} // namespace plates
} // namespace terraneo
/*
* Copyright (c) 2022 Berta Vilacis, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program 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.
*
* 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.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "core/extern/json.hpp"
#include "terraneo/helpers/conversions.hpp"
#include "terraneo/helpers/typeAliases.hpp"
#include "terraneo/plates/functionsForRotations.hpp"
#include "terraneo/plates/types.hpp"
namespace terraneo {
namespace plates {
/// Class for managing plate topology information
class PlateStorage
{
public:
/// Information describing a single plate at a certain age
struct PlateInfo
{
/// Set to true, after plate was rotated to xy-plane
bool rotatedToXY{ false };
/// ID of the plate
uint_t id{ 0 };
/// Nodes describing the boundary of the plate
Polygon boundary;
/// Center of the plate in original coordinates
///
/// If the plate was a 2D object, this would represent its barycenter.
/// The center is not rotated to the xy-plane, but kept at its original
/// position.
vec3D center{ real_c( 0 ), real_c( 0 ), real_c( 0 ) };
/// Textual name of plate
std::string name;
};
using plateVec_t = std::vector< PlateInfo >;
using ageToPlatesMap_t = std::map< std::string, plateVec_t >;
template < typename ImportStrategy >
PlateStorage( std::string nameOfPlateTopologiesFile, real_t sphereRadius, ImportStrategy readJsonfile )
: srcFile_( nameOfPlateTopologiesFile )
, sphereRadius_( sphereRadius )
{
// import topology data
nlohmann::json rootNode = readJsonfile( nameOfPlateTopologiesFile )["ages"];
// convert required data into our internal format
extractPlateInfo( rootNode );
// current approach treats plates as approximately being 2D
rotatePlatesToXYPlane();
// for testing
printStatistics();
}
/// Report statistics on what an object of this class stores
void printStatistics()
{
uint_t nPlates{ 0 };
for ( auto entry : ageToPlatesMap_ )
{
nPlates += entry.second.size();
}
std::cout << "PlateStorage object:\n"
<< " - stores plates for " << ageToPlatesMap_.size() << " age stages\n"
<< " - stores a total of " << nPlates << " plates\n"
<< " - data was obtained from file = '" << srcFile_ << "'" << std::endl;
}
plateVec_t& getPlatesForStage( real_t age )
{
auto iter = ageToPlatesMap_.find( ageToKey( age ) );
if ( iter == ageToPlatesMap_.end() )
{
std::cerr << "No plates found for " << ageToKey( age ) << std::endl;
std::abort();
}
return iter->second;
}
const plateVec_t& getPlatesForStage( real_t age ) const
{
auto iter = ageToPlatesMap_.find( ageToKey( age ) );
if ( iter == ageToPlatesMap_.end() )
{
std::cerr << "No plates found for " << ageToKey( age ) << std::endl;
std::abort();
}
return iter->second;
}
private:
// assemble key from age
inline std::string ageToKey( real_t age ) const
{
std::stringstream key;
key.precision( 4 );
key << std::fixed;
key << "topology_" << age << "Ma_polygon";
return key.str();
}
/// method for data preparation
///
/// this method will convert the imported data into a format more suitable for
/// using it in the computations; it will convert the plate polygons into
/// cartesian coordinates, compute the barycenter of each plate and group
/// plates into vectors depending on their age stage
void extractPlateInfo( const nlohmann::json& rootNode )
{
for ( uint_t idx = 0; idx < rootNode.size(); ++idx )
{
// prepare key for map entry and vector to hold plates for this age stage
std::string ageName = rootNode[idx]["name"];
uint_t nPlates = rootNode[idx]["features"].size();
ageToPlatesMap_.emplace( ageName, plateVec_t( nPlates ) );
plateVec_t& plates = ageToPlatesMap_[ageName];
// extract plates for this age stage and put into vector of plates
for ( uint_t k = 0; k < nPlates; ++k )
{
// get ID and name
plates[k].id = rootNode[idx]["features"][k]["properties"]["PLATEID1"];
plates[k].name = rootNode[idx]["features"][k]["properties"]["NAME"];
// extract coordinates of polygon (order: longitude, latitude)
const nlohmann::json coordinates = rootNode[idx]["features"][k]["geometry"]["coordinates"][0];
for ( auto& element : coordinates )
{
plates[k].boundary.push_back( { element[0], element[1], real_c( 0 ) } );
}
// convert to cartesian coordinates and compute center of plate
for ( auto& node : plates[k].boundary )
{
node = terraneo::conversions::sph2cart( { node[0], node[1] }, sphereRadius_ );
plates[k].center += node;
}
plates[k].center /= plates[k].boundary.size();
}
}
};
/// Rotate plates to xy-plane, so that their pseudo-barycenter lies at the origin
void rotatePlatesToXYPlane()
{
for ( auto& mapElem : ageToPlatesMap_ )
{
plateVec_t& plates = mapElem.second;
for ( auto& plate : plates )
{
mat3D rotMtx = terraneo::plates::getRotationMatrixPolygon( plate.center );
for ( auto& node : plate.boundary )
{
node = rotMtx * node;
}
plate.rotatedToXY = true;
}
}
};
/// name of datafile from which object obtained information
std::string srcFile_;
/// radius of sphere on which the plates are located
real_t sphereRadius_;
/// map to allow retrieving all plates of a certain age stage by giving that age
ageToPlatesMap_t ageToPlatesMap_;
};
} // namespace plates
} // namespace terraneo
/*
* Copyright (c) 2022 Berta Vilacis, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program 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.
*
* 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.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "terraneo/dataimport/io.hpp"
#include "terraneo/plates/PlateRotationProvider.hpp"
#include "terraneo/plates/PlateStorage.hpp"
#include "terraneo/plates/SmoothingStrategies.hpp"
// preserve ordering of includes
#include "terraneo/plates/functionsForPlates.hpp"
namespace terraneo {
namespace plates {
/// API class for computation of velocities from plate reconstructions
class PlateVelocityProvider
{
public:
PlateVelocityProvider( std::string nameOfTopologiesFile, std::string nameOfRotationsFile )
: plateTopologies_( nameOfTopologiesFile,
real_c( 1 ),
[]( const std::string& filename ) { return terraneo::io::readJsonFile( filename ); } )
, plateRotations_( nameOfRotationsFile,
[]( const std::string& filename ) { return terraneo::io::readRotationsFile( filename ); } )
{}
vec3D getPointVelocity( const vec3D& point, const real_t age )
{
return getPointVelocity( point, age, LinearDistanceSmoother{ 0.015 } );
}
template < typename SmoothingStrategy >
vec3D getPointVelocity( const vec3D& point, const real_t age, SmoothingStrategy computeSmoothing )
{
uint_t plateID{ 0 };
bool plateFound{ false };
real_t distance{ real_c( -1 ) };
std::tie( plateFound, plateID, distance ) = findPlateAndDistance( age, plateTopologies_, point );
if ( plateFound )
{
std::cout << "Point found on plate with ID = " << plateID << ", distance to boundary = " << distance << std::endl;
}
else
{
std::cout << "Point not found on any plate!" << std::endl;
}
real_t smoothingFactor = computeSmoothing( distance );
std::cout << "Smoothing Factor: " << smoothingFactor << "\n\n";
std::cout << "Plate ID: " << plateID << "\n\n";
if ( plateID == 0 )
{
std::cout << "No plate ID assigned (Plate ID =0). Velocity set to (0.0,0.0,0.0)" << std::endl;
// return{ real_c{0}, real_c{0}, real_c{0} };
return { 0.0, 0.0, 0.0 };
}
else
{
return computeCartesianVelocityVector( plateRotations_, plateID, age, point, smoothingFactor );
}