-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathExchangeHalos.cpp
150 lines (136 loc) · 7.16 KB
/
ExchangeHalos.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
// Example Includes
#include "ExchangeHalos.hpp"
// C++ includes
#include <iostream>
#include <string>
/// @brief Evaluate some closed-form Spherical Harmonic functions with an optional multiplier term
/// @param lon Longitude in lat-lon space
/// @param lat Latitude in lat-lon space
/// @param type Function type
/// @param multiplier Optional multiplier to scale value (default=1.0)
/// @return
static double evaluate_function( double lon, double lat, int type = 1, double multiplier = 1.0 )
{
switch( type )
{
case 1:
return ( 2.0 + std::pow( sin( 2.0 * lat ), 16.0 ) * cos( 16.0 * lon ) ) * multiplier;
default:
return ( 2.0 + cos( lon ) * cos( lon ) * cos( 2.0 * lat ) ) * multiplier;
}
}
moab::ErrorCode RuntimeContext::create_sv_tags( moab::Tag& tagScalar, moab::Tag& tagVector,
moab::Range& entities ) const
{
// Get element (centroid) coordinates so that we can evaluate some arbitrary data
std::vector< double > entCoords = compute_centroids( entities ); // [entities * [lon, lat]]
if( proc_id == 0 ) std::cout << "> Getting scalar tag handle " << scalar_tagname << "..." << std::endl;
double defSTagValue = -1.0;
bool createdTScalar = false;
// Get or create the scalar tag: default name = "scalar_variable"
// Type: double, Components: 1, Layout: Dense (all entities potentially), Default: -1.0
runchk( moab_interface->tag_get_handle( scalar_tagname.c_str(), 1, moab::MB_TYPE_DOUBLE, tagScalar,
moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &defSTagValue, &createdTScalar ),
"Retrieving scalar tag handle failed" );
// we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
assert( createdTScalar );
// set the data for scalar tag with an analytical Spherical Harmonic function
{
std::vector< double > tagValues( entities.size(), -1.0 );
std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
static int index = 0;
const int offset = index++ * 2;
return evaluate_function( entCoords[offset], entCoords[offset + 1] );
} );
// Set local scalar tag data for exchange
runchk( moab_interface->tag_set_data( tagScalar, entities, tagValues.data() ),
"Setting scalar tag data failed" );
}
if( proc_id == 0 ) std::cout << "> Getting vector tag handle " << vector_tagname << "..." << std::endl;
std::vector< double > defVTagValue( vector_length, -1.0 );
bool createdTVector = false;
// Get or create the scalar tag: default name = "vector_variable"
// Type: double, Components: vector_length, Layout: Dense (all entities potentially), Default: [-1.0,..]
runchk( moab_interface->tag_get_handle( vector_tagname.c_str(), vector_length, moab::MB_TYPE_DOUBLE, tagVector,
moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, defVTagValue.data(),
&createdTVector ),
"Retrieving vector tag handle failed" );
// we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
assert( createdTVector );
// set the data for vector tag with an analytical Spherical Harmonic function
// with an optional scaling for each component; just to make it look different :-)
{
const int veclength = vector_length;
std::vector< double > tagValues( entities.size() * veclength, -1.0 );
std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
static int index = 0;
const int offset = ( index++ / veclength ) * 2;
return evaluate_function( entCoords[offset], entCoords[offset + 1], 2, ( index % veclength + 1.0 ) );
} );
// Set local tag data for exchange
runchk( moab_interface->tag_set_data( tagVector, entities, tagValues.data() ),
"Setting vector tag data failed" );
}
return moab::MB_SUCCESS;
}
moab::ErrorCode RuntimeContext::load_file( bool load_ghosts )
{
/// Parallel Read options:
/// PARALLEL = type {READ_PART} : Read on all tasks
/// PARTITION_METHOD = RCBZOLTAN : Use Zoltan partitioner to compute an online partition and redistribute on the
/// fly PARTITION = PARALLEL_PARTITION : Partition as you read based on part information stored in h5m file
/// PARALLEL_RESOLVE_SHARED_ENTS : Communicate to all processors to get the shared adjacencies
/// consistently in parallel
/// PARALLEL_GHOSTS : a.b.c
/// : a = 2 - highest dimension of entities (2D in this case)
/// : b = 1 - dimension of entities to calculate adjacencies (vertex=0, edges=1)
/// : c = 3 - number of ghost layers needed (3 in this case)
std::string read_options = "DEBUG_IO=0;";
std::string::size_type idx = input_filename.rfind( '.' );
std::string extension = "";
if( num_procs > 1 && idx != std::string::npos )
{
extension = input_filename.substr( idx + 1 );
if( !extension.compare( "nc" ) )
// PARTITION_METHOD= [RCBZOLTAN, TRIVIAL]
read_options += "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
"PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
else if( !extension.compare( "h5m" ) )
read_options +=
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;"
"PARALLEL_RESOLVE_SHARED_ENTS;" +
( load_ghosts ? "PARALLEL_THIN_GHOST_LAYER;PARALLEL_GHOSTS=2.1." + std::to_string( ghost_layers ) + ";"
: "" );
else
{
std::cout << "Error unsupported file type (only h5m and nc) for this example: " << input_filename
<< std::endl;
return moab::MB_UNSUPPORTED_OPERATION;
}
}
// Load the file from disk with given read options in parallel and associate all entities to fileset
return moab_interface->load_file( input_filename.c_str(), &fileset, read_options.c_str() );
}
std::vector< double > RuntimeContext::compute_centroids( const moab::Range& entities ) const
{
double node[3];
std::vector< double > eCentroids( entities.size() * 2 ); // [lon, lat]
size_t offset = 0;
for( auto entity : entities )
{
// Get the element coordinates (centroid) on the real mesh
runchk_cont( moab_interface->get_coords( &entity, 1, node ), "Getting entity coordinates failed" );
// scale by magnitude so that element is on unit sphere
double magnitude = std::sqrt( node[0] * node[0] + node[1] * node[1] + node[2] * node[2] );
node[0] /= magnitude;
node[1] /= magnitude;
node[2] /= magnitude;
// compute the spherical transformation onto unit sphere
eCentroids[offset] = atan2( node[1], node[0] );
if( eCentroids[offset] < 0.0 ) eCentroids[offset] += 2.0 * M_PI;
eCentroids[offset + 1] = asin( node[2] );
offset += 2; // increment the offset
}
// return centroid list for elements
return eCentroids;
}