Merge pull request #161 from AetherModel/develop

Develop

include/aether.h

@@ -80,9 +80,11 @@ using json = nlohmann::json;
 #include "cubesphere.h"
 #include "sphere.h"
 #include "quadtree.h"
+#include "dipole.h"
 
 // not done
 #include "grid.h"
+#include "init_mag_grid.h"
 
 // Contains the neutral states and derived quantities
 #include "neutrals.h"

include/dipole.h

@@ -0,0 +1,129 @@
+// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members)
+// Full license can be found in License.md
+
+#ifndef INCLUDE_DIPOLE_H_
+#define INCLUDE_DIPOLE_H_
+
+#include "aether.h"
+#include <memory>
+
+/*************************************************
+ * \brief A namespace with all (1-root) dipole grid logic.
+ *************************************************/
+namespace Dipole {
+
+  /// The normalized origins of each face of the cube (i.e. corner)
+  static const arma_mat ORIGINS = {
+				   { 0.0, -0.5, 0.0}
+  };
+
+  /// Normalized right steps in cube
+  static const arma_mat RIGHTS = {
+				  {2.0, 0.0, 0.0}
+  };
+
+  /// Normalized up steps in cube
+  static const arma_mat UPS = {
+			       {0.0, 1.0, 0.0}
+  };
+
+};
+
+/*************************************************
+ * \brief A namespace with all (2-root) dipole grid logic.
+ *************************************************/
+namespace Dipole2 {
+
+  /// The normalized origins of each face of the cube (i.e. corner)
+  static const arma_mat ORIGINS = {
+				   { 0.0, -0.5, 0.0},
+           { 0.0,  0.0, 0.0}
+  };
+
+  /// Normalized right steps in cube
+  static const arma_mat RIGHTS = {
+				  {2.0, 0.0, 0.0},
+				  {2.0, 0.0, 0.0}
+  };
+
+  /// Normalized up steps in cube
+  static const arma_mat UPS = {
+			       {0.0, 0.5, 0.0},
+			       {0.0, 0.5, 0.0}
+  };
+
+};
+
+/*************************************************
+ * \brief A namespace with all (4-root) dipole grid logic.
+ *************************************************/
+namespace Dipole4 {
+
+  /// The normalized origins of each face of the cube (i.e. corner)
+  static const arma_mat ORIGINS = {
+				   { 0.0, -0.5, 0.0},
+				   { 1.0, -0.5, 0.0},
+				   { 1.0,  0.0, 0.0},
+           { 0.0,  0.0, 0.0}
+  };
+
+  /// Normalized right steps in cube
+  static const arma_mat RIGHTS = {
+				  {1.0, 0.0, 0.0},
+				  {1.0, 0.0, 0.0},
+				  {1.0, 0.0, 0.0},
+				  {1.0, 0.0, 0.0}
+  };
+
+  /// Normalized up steps in cube
+  static const arma_mat UPS = {
+			       {0.0, 0.5, 0.0},
+			       {0.0, 0.5, 0.0},
+			       {0.0, 0.5, 0.0},
+			       {0.0, 0.5, 0.0}
+  };
+
+};
+
+
+/*************************************************
+ * \brief A namespace with all (6-root) dipole grid logic.
+ * This is the same as the Sphere6
+ *************************************************/
+
+namespace Dipole6 {
+
+/// The normalized origins of each face of the cube (i.e. corner)
+static const arma_mat ORIGINS = {
+    {    0.0, -0.5, 0.0},
+    {2.0/3.0, -0.5, 0.0},
+    {4.0/3.0, -0.5, 0.0},
+    {    0.0, 0.0, 0.0},
+    {2.0/3.0, 0.0, 0.0},
+    {4.0/3.0, 0.0, 0.0}
+};
+
+/// Normalized right steps in cube
+static const arma_mat RIGHTS = {
+		{ 2.0/3.0,  0.0, 0.0},
+		{ 2.0/3.0,  0.0, 0.0},
+		{ 2.0/3.0,  0.0, 0.0},
+		{ 2.0/3.0,  0.0, 0.0},
+		{ 2.0/3.0,  0.0, 0.0},
+		{ 2.0/3.0,  0.0, 0.0}
+};
+
+/// Normalized right steps in cube
+static const arma_mat UPS = {
+		{ 0.0, 0.5, 0.0},
+		{ 0.0, 0.5, 0.0},
+		{ 0.0, 0.5, 0.0},
+		{ 0.0, 0.5, 0.0},
+		{ 0.0, 0.5, 0.0},
+		{ 0.0, 0.5, 0.0}
+};
+
+}
+
+
+#endif  // INCLUDE_DIPOLE_H_

include/grid.h

@@ -73,10 +73,37 @@ class Grid
   arma_cube magAlt_scgc, magZ_scgc;
   arma_cube magLocalTime_scgc;
 
+  // Dipole coordinates:
+  // Phi => Longitude
+  // P   => L-shell
+  // Q   => Distance along field line
   arma_cube magPhi_scgc;
   arma_cube magP_scgc;
   arma_cube magQ_scgc;
 
+  // And the corners/edges for the magnetic grid:
+  arma_cube magLon_Left;
+  arma_cube magLon_Corner;
+
+  arma_cube magLat_Down;
+  arma_cube magLat_Below;
+  arma_cube magLat_Corner;
+
+  arma_cube magAlt_Down;
+  arma_cube magAlt_Below;
+  arma_cube magAlt_Corner;
+
+  //For easier interpolation:
+  arma_vec baseLats_down;
+
+  // these need to be stored in (p,q) coords for a bit, its messy:
+  arma_cube magP_Down;
+  arma_cube magP_Below;
+  arma_cube magQ_Down;
+  arma_cube magQ_Below;
+  arma_cube magP_Corner;
+  arma_cube magQ_Corner;
+
   // These are the locations of the magnetic poles:
   //  ll -> lat, lon, radius independent
   arma_vec mag_pole_north_ll;
@@ -213,6 +240,7 @@ class Grid
   void calc_alt_grid_spacing();
   void calc_lat_grid_spacing();
   void calc_long_grid_spacing();
+
   void fill_grid_radius(Planets planet);
   void calc_rad_unit(Planets planet);
   void calc_gravity(Planets planet);
@@ -230,21 +258,33 @@ class Grid
   void calc_cent_acc(Planets planet);
 
   // Make mag-field grid:
-  std::pair<precision_t, precision_t> lshell_to_qn_qs(Planets planet,
-                                                      precision_t Lshell,
-                                                      precision_t Lon,
-                                                      precision_t AltMin);
   void convert_dipole_geo_xyz(Planets planet, precision_t XyzDipole[3],
                               precision_t XyzGeo[3]);
-  std::pair<arma_vec, arma_vec> fill_dipole_q_line(precision_t qN_,
-                                                   precision_t qS_,
-                                                   precision_t Gamma_,
-                                                   int64_t nZ_,
-                                                   precision_t Lshell_,
-                                                   precision_t min_alt_);
-  std::pair<precision_t, precision_t> qp_to_r_theta(precision_t q, precision_t p);
-  void init_dipole_grid(Quadtree quadtree, Planets planet);
-  arma_vec rNorm1d, lat1dalong;
+
+  bool init_dipole_grid(Quadtree quadtree_ion, Planets planet);
+  // Support functions:
+  void calc_dipole_grid_spacing(Planets planet);
+  void calc_alt_dipole_grid_spacing();
+  void calc_lat_dipole_grid_spacing();
+  void calc_long_dipole_grid_spacing();
+  void fill_field_lines(arma_vec baseLats, precision_t min_altRe,
+                        precision_t Gamma, Planets planet,
+                        bool isCorner);
+  void dipole_alt_edges(Planets planet, precision_t min_altRe);
+  // get the latitude spacing given the quadtree start & size, and the latitude limits
+  // extent: quadtree up
+  // origin: quadtree origin
+  // upper_lim: upper latitude limit (input)
+  // lower_lim: lower latitude limit (from min_apex)
+  // nLats: number of latitudes (nY)
+  // spacing_factor: (not supported yet), so always 1.0. Will adjust baselat spacing, eventually.
+  arma_vec baselat_spacing(precision_t extent,
+                          precision_t origin,
+                          precision_t upper_lim,
+                          precision_t lower_lim,
+                          // int16_t nLats,
+                          precision_t spacing_factor);
+
   // Update ghost cells with values from other processors
   void exchange(arma_cube &data, const bool pole_inverse);
 

include/init_mag_grid.h

@@ -0,0 +1,35 @@
+// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members)
+// Full license can be found in License.md
+
+
+#ifndef INCLUDE_INIT_MAG_GRID_H_
+#define INCLUDE_INIT_MAG_GRID_H_
+
+#include "aether.h"
+#include "planets.h"
+#include "grid.h"
+
+bool init_dipole_grid(Grid &mGrid, Planets planet);
+
+// Analytic solution to get from q,p dipole coords to r,theta
+// q coordinate along b-field line
+// p l-shell 
+// return (r,theta)
+std::pair<precision_t, precision_t> qp_to_r_theta(precision_t q, precision_t p);
+
+// Take limits & specs of field line, fill it with points.
+// std:: tuple <arma_mat, arma_mat, arma_vec> Grid::fill_field_lines (arma_vec lShells, 
+//                                                 int64_t nAlts,
+//                                                 int64_t nLats,
+//                                                 precision_t Gamma);
+
+
+// convert mag to geographic
+std::vector <arma_cube> mag_to_geo(arma_cube magLon, 
+                                   arma_cube magLat,
+                                   arma_cube magAlt,
+                                   Planets planet);
+
+
+#endif  // INCLUDE_INIT_GEO_GRID_H_
+

include/inputs.h

@@ -63,7 +63,7 @@ class Inputs {
     // Only needed for Mag Field grid:
     // min_apex (not used) and LatStretch is used 
     // as lat = min_lat + dlat where dlat = acos(cos(lat^stretch))^(1/stretch)
-    precision_t min_apex, LatStretch, FieldLineStretch, max_lat_dipole;
+    precision_t min_apex, LatStretch, FieldLineStretch, max_blat;
 
     // Some grid shapes allow specification of altitudes based on a file:
     std::string alt_file;

include/parallel.h

@@ -34,7 +34,7 @@ extern MPI_Comm aether_comm;
 /**********************************************************************
   \brief initialize mpi and figure out ensembles and grid blocks
 **/
-bool init_parallel(Quadtree &quadtree);
+bool init_parallel(Quadtree &quadtree, Quadtree &quadtree_ion);
 
 /**********************************************************************
   \brief Pack variables for message passing

include/sphere.h

@@ -41,7 +41,7 @@ static const arma_mat ORIGINS = {
     {4.0/3.0, -0.5, 0.0},
     {    0.0, 0.0, 0.0},
     {2.0/3.0, 0.0, 0.0},
-    {3.0/3.0, 0.0, 0.0}
+    {4.0/3.0, 0.0, 0.0}
 };
 
 /// Normalized right steps in cube

share/run/UA/inputs/defaults.json

@@ -80,22 +80,16 @@
         "AltFile" : ""},
 
     "ionGrid" : {
-        "Shape" : "sphere",
-	      "LatRange" : [-90.0, 90.0],
+        "Shape" : "dipole",
         "nLatsPerBlock" : 18,
-	      "LonRange" : [0.0, 360.0],
+		"LonRange" : [0.0, 360.0],
         "nLonsPerBlock" : 22,
         "nAlts" : 50,
         "MinAlt" : 80.0,
         "MinApex" : 120.0,
-        "MaxAlt" : 5000.0,
         "LatMax":88.0,
-        "LineSpacing":0.2,
-        "LatStretch":6.0,
-        "dAltkm" : 5.0,
-        "dAltScale" : 0.25,
-	      "IsUniformAlt" : true,
-        "AltFile" : ""},
+        "LatStretch":1.0,
+        "dAltStretch" : 0.6},
 
     "Oblate" : {
 	"isOblate" : false,

share/run/aether.json

@@ -26,11 +26,12 @@
         "IsUniformAlt" : false},
 
     "ionGrid": {
-        "dAltScale": 0.2,
-        "LatStretch": 6,
+        "dAltStretch": 0.6,
+        "LatStretch": 1,
         "Shape": "dipole",
-        "nLatsPerBlock" : 22,
-        "nAlts":159,
+        "nLonsPerBlock": 36,
+        "nLatsPerBlock" : 18,
+        "nAlts":36,
         "LatMax":88,
         "MinAlt":100.0,
         "MinApex":150.0