forked from fmihpc/vlasiator
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathobject_wrapper.cpp
215 lines (179 loc) · 12 KB
/
object_wrapper.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
#include "object_wrapper.h"
#include "velocity_mesh_parameters.h"
#include "readparameters.h"
#include <iostream>
bool ObjectWrapper::addParameters() {
typedef Readparameters RP;
// Parameters needed to create particle populations
if (RP::helpRequested) { // dummy name for the help message
RP::add("ParticlePopulations","Name of the simulated particle populations (string)","<population>");
} else {
RP::addComposing("ParticlePopulations","Name of the simulated particle populations (string)");
}
return true;
}
bool ObjectWrapper::addPopulationParameters() {
typedef Readparameters RP;
std::vector<std::string> popNames;
if (RP::helpRequested) {
popNames.push_back(std::string("<population>"));
} else {
RP::get("ParticlePopulations", popNames);
}
// Create appropriate subparameters for each population
for(auto& pop : popNames) {
species::Species newSpecies;
vmesh::MeshParameters newVMesh;
// Originally, there was support for species and velocity meshes to be separate.
// This was abandoned, since there wasn't really any use for it.
newSpecies.name = newVMesh.name = pop;
newSpecies.velocityMesh = vmesh::getMeshWrapper()->velocityMeshesCreation->size();
getObjectWrapper().particleSpecies.push_back(newSpecies);
vmesh::getMeshWrapper()->velocityMeshesCreation->push_back(newVMesh);
RP::add(pop + "_properties.charge", "Particle charge, in units of elementary charges (int)", 1);
RP::add(pop + "_properties.mass_units", "Units in which particle mass is given, either 'PROTON' or 'ELECTRON' (string)", std::string("PROTON"));
RP::add(pop + "_properties.mass","Particle mass in given units (float)", 1);
// Grid sparsity parameters
RP::add(pop + "_sparse.minValue", "Minimum value of distribution function in any cell of a velocity block for the block to be considered to have contents", 1e-15);
RP::add(pop + "_sparse.blockAddWidthV", "Number of layers of blocks that are kept in velocity space around the blocks with content",1);
RP::add(pop + "_sparse.conserve_mass", "If true, then mass is conserved by scaling the dist. func. in the remaining blocks", false);
RP::add(pop + "_sparse.dynamicAlgorithm", "Type of algorithm used for calculating the dynamic minValue; 0 = none, 1 = linear algorithm based on rho, 2 = linear algorithm based on Blocks, (Example linear algorithm: y = kx+b, where dynamicMinValue1=k*dynamicBulkValue1 + b, and dynamicMinValue2 = k*dynamicBulkValue2 + b", 0);
RP::add(pop + "_sparse.dynamicMinValue1", "The minimum value for the dynamic minValue", 1);
RP::add(pop + "_sparse.dynamicMinValue2", "The maximum value (value 2) for the dynamic minValue", 1);
RP::add(pop + "_sparse.dynamicBulkValue1", "Minimum value for the dynamic algorithm range, so for example if dynamicAlgorithm=1 then for sparse.dynamicBulkValue1 = 1e3, sparse.dynamicBulkValue2=1e5, we apply the algorithm to cells for which 1e3<cell.rho<1e5", 0);
RP::add(pop + "_sparse.dynamicBulkValue2", "Maximum value for the dynamic algorithm range, so for example if dynamicAlgorithm=1 then for sparse.dynamicBulkValue1 = 1e3, sparse.dynamicBulkValue2=1e5, we apply the algorithm to cells for which 1e3<cell.rho<1e5", 0);
// Grid parameters
RP::add(pop + "_vspace.vx_min","Minimum value for velocity mesh vx-coordinates.",0);
RP::add(pop + "_vspace.vx_max","Maximum value for velocity mesh vx-coordinates.",0);
RP::add(pop + "_vspace.vy_min","Minimum value for velocity mesh vy-coordinates.",0);
RP::add(pop + "_vspace.vy_max","Maximum value for velocity mesh vx-coordinates.",0);
RP::add(pop + "_vspace.vz_min","Minimum value for velocity mesh vz-coordinates.",0);
RP::add(pop + "_vspace.vz_max","Maximum value for velocity mesh vx-coordinates.",0);
RP::add(pop + "_vspace.vx_length","Initial number of velocity blocks in vx-direction.",1);
RP::add(pop + "_vspace.vy_length","Initial number of velocity blocks in vy-direction.",1);
RP::add(pop + "_vspace.vz_length","Initial number of velocity blocks in vz-direction.",1);
RP::add(pop + "_vspace.max_refinement_level","Maximum allowed mesh refinement level.", 1);
// Thermal / suprathermal parameters
Readparameters::add(pop + "_thermal.vx", "Center coordinate for the maxwellian distribution. Used for calculating the suprathermal moments.", -500000.0);
Readparameters::add(pop + "_thermal.vy", "Center coordinate for the maxwellian distribution. Used for calculating the suprathermal moments.", 0.0);
Readparameters::add(pop + "_thermal.vz", "Center coordinate for the maxwellian distribution. Used for calculating the suprathermal moments.", 0.0);
Readparameters::add(pop + "_thermal.radius", "Radius of the maxwellian distribution. Used for calculating the suprathermal moments. If set to 0 (default), the thermal/suprathermal DROs are skipped.", 0.0);
// Precipitation parameters
Readparameters::add(pop + "_precipitation.nChannels", "Number of energy channels for precipitation differential flux evaluation", 16);
Readparameters::add(pop + "_precipitation.emin", "Lowest energy channel (in eV) for precipitation differential flux evaluation", 0.1);
Readparameters::add(pop + "_precipitation.emax", "Highest energy channel (in eV) for precipitation differential flux evaluation", 100.0);
Readparameters::add(pop + "_precipitation.lossConeAngle", "Fixed loss cone opening angle (in deg) for precipitation differential flux evaluation", 10.0);
// Energy density parameters
Readparameters::add(pop + "_energydensity.limit1", "Lower limit of second bin for energy density, given in units of solar wind ram energy.", 5.0);
Readparameters::add(pop + "_energydensity.limit2", "Lower limit of third bin for energy density, given in units of solar wind ram energy.", 10.0);
Readparameters::add(pop + "_energydensity.solarwindspeed", "Incoming solar wind velocity magnitude in m/s. Used for calculating energy densities.", 0.0);
Readparameters::add(pop + "_energydensity.solarwindenergy", "Incoming solar wind ram energy in eV. Used for calculating energy densities.", 0.0);
}
return true;
}
bool ObjectWrapper::getPopulationParameters() {
typedef Readparameters RP;
// Particle population parameters
for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) {
species::Species& species=getObjectWrapper().particleSpecies[i];
vmesh::MeshParameters& vMesh=vmesh::getMeshWrapper()->velocityMeshesCreation->at(i);
const std::string& pop = species.name;
// Sanity check name
if(species.name != vMesh.name) {
std::cerr << "ParticlePopulation parse error: Name " << species.name << " != " << vMesh.name << std::endl;
return false;
}
// Elementary particle parameters
RP::get(pop + "_properties.charge", species.charge);
species.charge *= physicalconstants::CHARGE;
RP::get(pop + "_properties.mass", species.mass);
std::string massUnit;
RP::get(pop + "_properties.mass_units", massUnit);
if(massUnit == "PROTON") {
species.mass *= physicalconstants::MASS_PROTON;
} else if(massUnit == "ELECTRON") {
species.mass *= physicalconstants::MASS_ELECTRON;
} else {
std::cerr << "Invalid mass unit for species " << pop << ": '" << massUnit << "'" << std::endl;
return false;
}
// sparsity parameters
RP::get(pop + "_sparse.minValue", species.sparseMinValue);
RP::get(pop + "_sparse.blockAddWidthV", species.sparseBlockAddWidthV);
RP::get(pop + "_sparse.conserve_mass", species.sparse_conserve_mass);
RP::get(pop + "_sparse.dynamicAlgorithm", species.sparseDynamicAlgorithm);
RP::get(pop + "_sparse.dynamicBulkValue1", species.sparseDynamicBulkValue1);
RP::get(pop + "_sparse.dynamicBulkValue2", species.sparseDynamicBulkValue2);
RP::get(pop + "_sparse.dynamicMinValue1", species.sparseDynamicMinValue1);
RP::get(pop + "_sparse.dynamicMinValue2", species.sparseDynamicMinValue2);
// Particle velocity space properties
RP::get(pop + "_vspace.vx_min",vMesh.meshLimits[0]);
RP::get(pop + "_vspace.vx_max",vMesh.meshLimits[1]);
RP::get(pop + "_vspace.vy_min",vMesh.meshLimits[2]);
RP::get(pop + "_vspace.vy_max",vMesh.meshLimits[3]);
RP::get(pop + "_vspace.vz_min",vMesh.meshLimits[4]);
RP::get(pop + "_vspace.vz_max",vMesh.meshLimits[5]);
RP::get(pop + "_vspace.vx_length",vMesh.gridLength[0]);
RP::get(pop + "_vspace.vy_length",vMesh.gridLength[1]);
RP::get(pop + "_vspace.vz_length",vMesh.gridLength[2]);
if(vMesh.gridLength[0] > MAX_BLOCKS_PER_DIM ||
vMesh.gridLength[1] > MAX_BLOCKS_PER_DIM ||
vMesh.gridLength[2] > MAX_BLOCKS_PER_DIM ) {
// Build error message as a string first, so that the cerr output hapens atomically and we don't spam
// thousands of unreadable lines through each other
std::string errormsg = "(VSPACE) ERROR: Velocity mesh for population " + species.name + " has too many blocks per dimension. Maximum defined in MAX_BLOCKS_PER_DIM is " + std::to_string(MAX_BLOCKS_PER_DIM) + " "
+ std::string(__FILE__) + ":" + std::to_string(__LINE__) + "\n";
std::cerr << errormsg;
}
/* Special handling of GPU version with WID=8; halve the number of blocks */
#ifdef USE_GPU
int myRank;
MPI_Comm_rank(MPI_COMM_WORLD,&myRank);
if ((WID==8 && P::adaptGPUWID)) {
// First verify that we can halve the value2
if ( (vMesh.gridLength[0]%2==0) && (vMesh.gridLength[1]%2==0) && (vMesh.gridLength[2]%2==0)) {
vMesh.gridLength[0] /= 2;
vMesh.gridLength[1] /= 2;
vMesh.gridLength[2] /= 2;
if(myRank==MASTER_RANK) {
std::cerr<<" Note: Using GPU mode with WID=8; Halving velocity block counts per dimension. Deactivate with parameter adaptGPUWID=false."<<std::endl;
}
} else {
if(myRank==MASTER_RANK) {
std::cerr<<" Warning: Using GPU mode with WID=8 but odd number of velocity blocks! Cannot halve the blocks count."<<std::endl;
}
}
}
#endif
vMesh.blockLength[0] = vMesh.blockLength[1] = vMesh.blockLength[2] = WID;
int maxRefLevel; // Temporary variable, since target value is a uint8_t
RP::get(pop + "_vspace.max_refinement_level",maxRefLevel);
vMesh.refLevelMaxAllowed = maxRefLevel;
//Get thermal / suprathermal moments parameters
Readparameters::get(pop + "_thermal.radius", species.thermalRadius);
Readparameters::get(pop + "_thermal.vx", species.thermalV[0]);
Readparameters::get(pop + "_thermal.vy", species.thermalV[1]);
Readparameters::get(pop + "_thermal.vz", species.thermalV[2]);
//Get energy density parameters
Readparameters::get(pop + "_energydensity.limit1", species.EnergyDensityLimit1);
Readparameters::get(pop + "_energydensity.limit2", species.EnergyDensityLimit2);
Readparameters::get(pop + "_energydensity.solarwindenergy", species.SolarWindEnergy);
Readparameters::get(pop + "_energydensity.solarwindspeed", species.SolarWindSpeed);
const Real EPSILON = 1.e-25;
if (species.SolarWindEnergy < EPSILON) {
// Energy stored internally in SI units
species.SolarWindEnergy = 0.5 * species.mass * species.SolarWindSpeed * species.SolarWindSpeed;
} else {
species.SolarWindEnergy = species.SolarWindEnergy*physicalconstants::CHARGE;
}
// Get precipitation parameters
Readparameters::get(pop + "_precipitation.nChannels", species.precipitationNChannels);
Readparameters::get(pop + "_precipitation.emin", species.precipitationEmin);
Readparameters::get(pop + "_precipitation.emax", species.precipitationEmax);
Readparameters::get(pop + "_precipitation.lossConeAngle", species.precipitationLossConeAngle);
// Convert from eV to SI units
species.precipitationEmin = species.precipitationEmin*physicalconstants::CHARGE;
species.precipitationEmax = species.precipitationEmax*physicalconstants::CHARGE;
}
return true;
}