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

Makes PlateVelocityDemo more flexible and talkative

parent 259c9822
Pipeline #41047 passed with stages
in 121 minutes and 6 seconds
......@@ -42,32 +42,67 @@ using walberla::real_t;
using namespace hyteg;
using namespace terraneo;
typedef enum
{
PLATE_IDS,
VELOCITIES,
VELOCITIES_AND_IDS
} job_t;
std::string genOutputFileName( real_t age )
{
std::stringstream sstr;
sstr << "PlateVelocitiesDemo_age=";
sstr << std::setfill( '0' );
sstr << std::setw( 5 );
sstr << std::fixed;
sstr << std::setprecision( 1 );
sstr << age;
return sstr.str();
}
std::tuple< uint_t, uint_t > decodeRangeString( std::string& rangeStr )
{
size_t fPos = rangeStr.find( '-' );
if ( fPos == std::string::npos )
{
WALBERLA_ABORT( "Cannot decode range string '" << rangeStr << "'" );
}
uint_t stageInit = stoul( rangeStr.substr( 0, fPos ) );
uint_t stageStop = stoul( rangeStr.substr( fPos + 1 ) );
return std::make_tuple( stageInit, stageStop );
}
// =================================================================================
// Function to test computation of a velocity field from plate reconstruction data
// for all DoFs on the surface of a sphere
// =================================================================================
template < typename feFuncType >
void exportDemoFunc( uint_t level, std::shared_ptr< hyteg::PrimitiveStorage > storage, real_t age )
void performComputations( uint_t level,
std::shared_ptr< hyteg::PrimitiveStorage > storage,
real_t age,
job_t jobType,
terraneo::plates::PlateVelocityProvider& oracle,
std::string indent )
{
// initialise an oracle
std::string dataDir{"../../../data/terraneo/plates/"};
std::string fnameTopologies = dataDir + "topologies0-100Ma.geojson";
std::string fnameReconstructions = dataDir + "Global_EarthByte_230-0Ma_GK07_AREPS.rot";
terraneo::plates::PlateVelocityProvider oracle( fnameTopologies, fnameReconstructions );
// need that here, to capture it below ;-)
uint_t coordIdx = 0;
uint_t coordIdx = 0;
terraneo::plates::StatisticsPlateNotFoundHandler handlerWithStatistics;
std::function< real_t( const Point3D& ) > computeVelocityComponent = [&oracle, age, &coordIdx, &handlerWithStatistics]( const Point3D& point ) {
vec3D coords{point[0], point[1], point[2]};
vec3D velocity = oracle.getPointVelocity( coords, age, terraneo::plates::LinearDistanceSmoother{0.015}, handlerWithStatistics );
return velocity[int_c( coordIdx )];
};
// callback function for computing the velocity components
std::function< real_t( const Point3D& ) > computeVelocityComponent =
[&oracle, age, &coordIdx, &handlerWithStatistics]( const Point3D& point ) {
vec3D coords{ point[0], point[1], point[2] };
vec3D velocity =
oracle.getPointVelocity( coords, age, terraneo::plates::LinearDistanceSmoother{ 0.015 }, handlerWithStatistics );
return velocity[int_c( coordIdx )];
};
// callback function for determining plate IDs
std::function< real_t( const Point3D& ) > findPlateID = [&oracle, age]( const Point3D& point ) {
vec3D coords{point[0], point[1], point[2]};
vec3D coords{ point[0], point[1], point[2] };
return oracle.findPlateID( coords, age );
};
......@@ -76,26 +111,50 @@ void exportDemoFunc( uint_t level, std::shared_ptr< hyteg::PrimitiveStorage > st
BoundaryUID regionSurface = bcs.createDirichletBC( "surface", MeshInfo::flagOuterBoundary );
BoundaryUID regionCMB = bcs.createDirichletBC( "core-mantle-boundary", MeshInfo::flagInnerBoundary );
// set everything to zero in the interior and on the CMB
feFuncType demoFunc( "demoFunc", storage, level, level );
demoFunc.setBoundaryCondition( bcs );
demoFunc.interpolate( {real_c( 0 ), real_c( 0 ), real_c( 0 )}, level, Inner );
demoFunc.interpolate( {real_c( 0 ), real_c( 0 ), real_c( 0 )}, level, regionCMB );
// use pointers so that we only request memory for function in the desired jobType
std::shared_ptr< feFuncType > surfaceVelocity{ nullptr };
std::shared_ptr< typename feFuncType::VectorComponentType > plateID{ nullptr };
typename feFuncType::VectorComponentType plates( "plateID", storage, level, level );
plates.setBoundaryCondition( bcs );
plates.interpolate( real_c( 0 ), level, Inner );
plates.interpolate( real_c( 0 ), level, regionCMB );
// for DoF locations on the surface determine their associate plate IDs
if ( jobType == PLATE_IDS || jobType == VELOCITIES_AND_IDS )
{
plateID = std::make_shared< typename feFuncType::VectorComponentType >( "plateID", storage, level, level );
plateID->setBoundaryCondition( bcs );
plateID->interpolate( real_c( 0 ), level, Inner );
plateID->interpolate( real_c( 0 ), level, regionCMB );
plateID->interpolate( findPlateID, level, regionSurface );
}
for ( coordIdx = 0; coordIdx < 3; ++coordIdx )
// for DoF locations on the surface determine their associate velocity vectors
if ( jobType == VELOCITIES || jobType == VELOCITIES_AND_IDS )
{
demoFunc[coordIdx].interpolate( computeVelocityComponent, level, regionSurface );
surfaceVelocity = std::make_shared< feFuncType >( "plateVelocities", storage, level, level );
surfaceVelocity->setBoundaryCondition( bcs );
surfaceVelocity->interpolate( { real_c( 0 ), real_c( 0 ), real_c( 0 ) }, level, Inner );
surfaceVelocity->interpolate( { real_c( 0 ), real_c( 0 ), real_c( 0 ) }, level, regionCMB );
for ( coordIdx = 0; coordIdx < 3; ++coordIdx )
{
( *surfaceVelocity )[coordIdx].interpolate( computeVelocityComponent, level, regionSurface );
}
}
plates.interpolate( findPlateID, level, regionSurface );
hyteg::VTKOutput vtkOutput( "./output", "PlateVelocities", storage );
vtkOutput.add( plates );
vtkOutput.add( demoFunc );
// export results
std::string fName = genOutputFileName( age );
WALBERLA_LOG_INFO_ON_ROOT( "" << indent << "Exporting simulation data to file with basename '" << fName << "'" );
hyteg::VTKOutput vtkOutput( "./output", fName, storage );
switch ( jobType )
{
case PLATE_IDS:
vtkOutput.add( *plateID );
break;
case VELOCITIES:
vtkOutput.add( *surfaceVelocity );
break;
case VELOCITIES_AND_IDS:
vtkOutput.add( *plateID );
vtkOutput.add( *surfaceVelocity );
break;
}
vtkOutput.write( level, 0 );
handlerWithStatistics.generateReport();
......@@ -113,6 +172,7 @@ int main( int argc, char* argv[] )
// ------------
// Parameters
// ------------
WALBERLA_LOG_INFO_ON_ROOT( "*** STEP 1: Obtaining Steering Parameters" );
// check if a config was given on command line or load default file otherwise
auto cfg = std::make_shared< walberla::config::Config >();
......@@ -129,9 +189,40 @@ int main( int argc, char* argv[] )
const walberla::Config::BlockHandle params = cfg->getBlock( "Parameters" );
if ( walberla::MPIManager::instance()->worldRank() == 0 )
{
WALBERLA_LOG_INFO( "Running with the following steering parameters:" );
params.listParameters();
}
// determine job type
job_t jobType;
const std::string str = params.getParameter< std::string >( "jobType" );
if ( str == "plateIDs" )
{
jobType = PLATE_IDS;
}
else if ( str == "velocities" )
{
jobType = VELOCITIES;
}
else if ( str == "both" )
{
jobType = VELOCITIES_AND_IDS;
}
else
{
WALBERLA_ABORT( "Invalid jobType in parameter file!" );
}
// determine age indices
std::string rangeStr = params.getParameter< std::string >( "ageRange" );
auto [stageInit, stageStop] = decodeRangeString( rangeStr );
// ---------
// Meshing
// ---------
WALBERLA_LOG_INFO_ON_ROOT( "*** STEP 2: Generating Mesh" );
const uint_t level = params.getParameter< uint_t >( "level" );
const uint_t nRad = params.getParameter< uint_t >( "nRad" );
const uint_t nTan = params.getParameter< uint_t >( "nTan" );
......@@ -145,19 +236,44 @@ int main( int argc, char* argv[] )
std::shared_ptr< walberla::WcTimingTree > timingTree( new walberla::WcTimingTree() );
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage, timingTree );
// ============
// Delegation
// ============
std::string feSpace = params.getParameter< std::string >( "feSpace" );
const real_t age = params.getParameter< real_t >( "age" );
WALBERLA_LOG_INFO_ON_ROOT( "" << setupStorage );
// --------
// Oracle
// --------
WALBERLA_LOG_INFO_ON_ROOT( "*** STEP 3: Generating an Oracle" );
std::string dataDir{ "../../../data/terraneo/plates/" };
std::string fnameTopologies = dataDir + "topologies0-100Ma.geojson";
std::string fnameReconstructions = dataDir + "Global_EarthByte_230-0Ma_GK07_AREPS.rot";
terraneo::plates::PlateVelocityProvider oracle( fnameTopologies, fnameReconstructions );
if ( feSpace == "P1" )
WALBERLA_LOG_INFO_ON_ROOT( "*** STEP 4: Checking plate stages to work with" );
auto stages = oracle.getListOfPlateStages();
if ( stageStop >= stages.size() )
{
exportDemoFunc< P1VectorFunction< real_t > >( level, storage, age );
WALBERLA_ABORT( "ageRange parameter does not work for available plate stages!" );
}
else if ( feSpace == "P2" )
WALBERLA_LOG_INFO_ON_ROOT( " - Running from stage[" << stageInit << "] = " << stages[stageInit] << " Ma back to stage["
<< stageStop << "] = " << stages[stageStop] << " Ma" );
// ------------
// Delegation
// ------------
WALBERLA_LOG_INFO_ON_ROOT( "*** STEP 5: Running the actual computations" );
std::string feSpace = params.getParameter< std::string >( "feSpace" );
for ( uint_t currentStage = stageInit; currentStage <= stageStop; ++currentStage )
{
exportDemoFunc< P2VectorFunction< real_t > >( level, storage, age );
WALBERLA_LOG_INFO_ON_ROOT( " - stage[" << std::setw( 3 ) << currentStage << "] = " << std::fixed << std::setprecision( 1 )
<< stages[currentStage] << " Ma" );
if ( feSpace == "P1" )
{
performComputations< P1VectorFunction< real_t > >( level, storage, stages[currentStage], jobType, oracle, " " );
}
else if ( feSpace == "P2" )
{
performComputations< P2VectorFunction< real_t > >( level, storage, stages[currentStage], jobType, oracle, " " );
}
}
return EXIT_SUCCESS;
......
Parameters
{
// Decide what should be computed:
// plateIDs ...... determine plate IDs for surface points/DoFs
// velocities .... compute velocity vectors for surface points/DoFs
// both .......... do both
jobType plateIDs;
// The topology file 'topologies0-100Ma.geojson' the demo is using
// stores plates for 1001 age stages. You can specify a range of
// ages here for which you want to compute plate IDs and/or velocities.
// The range must be ints with 0 = 0 Ma and 1000 = 100 Ma
// (whitespace around - does not matter)
ageRange 27-30;
// meshing parameters
nTan 3;
nRad 2;
......@@ -7,7 +20,4 @@ Parameters
// FE space
feSpace P2;
// age stage
age 1.0;
}
......@@ -46,9 +46,9 @@ class PlateVelocityProvider
uint_t findPlateID( const vec3D& point, const real_t age )
{
uint_t plateID{idWhenNoPlateFound};
bool plateFound{false};
real_t distance{real_c( -1 )};
uint_t plateID{ idWhenNoPlateFound };
bool plateFound{ false };
real_t distance{ real_c( -1 ) };
std::tie( plateFound, plateID, distance ) = findPlateAndDistance( age, plateTopologies_, point );
......@@ -57,7 +57,7 @@ class PlateVelocityProvider
vec3D getPointVelocity( const vec3D& point, const real_t age )
{
return getPointVelocity( point, age, LinearDistanceSmoother{0.015}, DefaultPlateNotFoundHandler{} );
return getPointVelocity( point, age, LinearDistanceSmoother{ 0.015 }, DefaultPlateNotFoundHandler{} );
}
template < typename SmoothingStrategy, typename PlateNoteFoundStrategy >
......@@ -66,9 +66,9 @@ class PlateVelocityProvider
SmoothingStrategy computeSmoothing,
PlateNoteFoundStrategy&& errorHandler )
{
uint_t plateID{0};
bool plateFound{false};
real_t distance{real_c( -1 )};
uint_t plateID{ 0 };
bool plateFound{ false };
real_t distance{ real_c( -1 ) };
std::tie( plateFound, plateID, distance ) = findPlateAndDistance( age, plateTopologies_, point );
......@@ -88,7 +88,9 @@ class PlateVelocityProvider
return computeCartesianVelocityVector( plateRotations_, plateID, age, point, smoothingFactor );
};
const uint_t idWhenNoPlateFound{0};
const std::vector< real_t >& getListOfPlateStages() const { return plateTopologies_.getListOfPlateStages(); }
const uint_t idWhenNoPlateFound{ 0 };
private:
PlateStorage plateTopologies_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment