ionosphere.h

Go to the documentation of this file.

/*
 * This file is part of Vlasiator.
 * Copyright 2010-2020 Finnish Meteorological Institute
 *
 * For details of usage, see the COPYING file and read the "Rules of the Road"
 * at http://www.physics.helsinki.fi/vlasiator/
 *
 * 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 2 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, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */
 
#ifndef IONOSPHERE_H
#define IONOSPHERE_H
 
#include <cstdint>
#include <vector>
#include <functional>
#include "../definitions.h"
#include "../readparameters.h"
#include "../spatial_cells/spatial_cell_wrapper.hpp"
#include "sysboundarycondition.h"
#include "../backgroundfield/fieldfunction.hpp"
#include "../fieldsolver/fs_common.h"
 
#include <Eigen/Sparse>
#include <Eigen/Geometry>
 
using namespace projects;
using namespace std;
 
namespace SBC {
 
   // Hardcoded constants for calculating ion production table
   // TODO: Make these parameters?
   constexpr static int productionNumAccEnergies = 60;
   constexpr static int productionNumTemperatures = 60;
   constexpr static int productionNumParticleEnergies = 100;
   constexpr static Real productionMinAccEnergy = 0.1; // keV
   constexpr static Real productionMaxAccEnergy = 100.; // keV
   constexpr static Real productionMinTemperature = 0.1; // keV
   constexpr static Real productionMaxTemperature = 100.; // keV
   constexpr static Real ion_electron_T_ratio = 4.; // TODO: Make this a parameter (and/or find value from kinetics)
 
 
   struct IonosphereSpeciesParameters {
      Real rho;
      Real V0[3];
      Real T;
   };
 
   enum IonosphereBoundaryVDFmode { // How are inner boundary VDFs constructed from the ionosphere
      FixedMoments,      // Predefine temperature, density and V = EXB drift on the inner boundary.
      AverageMoments,    // Copy averaged density and temperature from nearest cells, V = EXB drift
      AverageAllMoments, // Same as above, but also copy V + add EXB drift to it
      CopyAndLosscone
   };
   extern IonosphereBoundaryVDFmode boundaryVDFmode;
 
   static const int MAX_TOUCHING_ELEMENTS = 12; // Maximum number of elements touching one node
   static const int MAX_DEPENDING_NODES = 22;   // Maximum number of depending nodes
 
   typedef Real iSolverReal; // Datatype for the ionosphere solver internal state
 
   // Ionosphere finite element grid
   struct SphericalTriGrid {
 
      // One finite element, spanned between 3 nodes
      struct Element {
         int refLevel = 0;
         std::array<uint32_t, 3> corners;                 // Node indices in the corners of this element
 
      };
      std::vector<Element> elements;
      std::vector<Eigen::Vector3d> elementCurlFreeCurrent;
      std::vector<Eigen::Vector3d> elementDivFreeCurrent;
      // One grid node
      struct Node {
         // Elements touching this node
         uint numTouchingElements=0;
         std::array<uint32_t, MAX_TOUCHING_ELEMENTS> touchingElements;
 
         // List of nodes the current node depends on
         uint numDepNodes = 0;
         std::array<uint32_t, MAX_DEPENDING_NODES> dependingNodes;
         std::array<Real, MAX_DEPENDING_NODES> dependingCoeffs;// Dependency coefficients
         std::array<Real, MAX_DEPENDING_NODES> transposedCoeffs; // Transposed dependency coefficient
 
         std::array<Real, 3> x = {0,0,0}; // Coordinates of the node
         std::array<Real, 3> xMapped = {0,0,0}; // Coordinates mapped along fieldlines into simulation domain
         int haveCouplingData = 0; // Does this rank carry coupling coordinate data for this node? (0 or 1)
         std::array<iSolverReal, N_IONOSPHERE_PARAMETERS> parameters = {0}; // Parameters carried by the node, see common.h
 
         int openFieldLine; 
 
         // Some calculation helpers
         Real electronDensity() { // Electron Density
            return parameters[ionosphereParameters::RHON];
         }
         Real electronTemperature() { // Electron Temperature
            return parameters[ionosphereParameters::TEMPERATURE];
         }
         Real deltaPhi() { // Field aligned potential drop between i'sphere and m'sphere
 
            // When the Knight-parameter is irrelevant, we can set this to zero
            return 0;
 
            // Alternative: Calculate it just like GUMCS does
 
            //if(electronDensity() == 0) {
            //   return 0;
            //}
 
            //Real retval = physicalconstants::K_B * electronTemperature() / physicalconstants::CHARGE
            //   * ((parameters[ionosphereParameters::SOURCE] / (physicalconstants::CHARGE * electronDensity()))
            //   * sqrt(2. * M_PI * physicalconstants::MASS_ELECTRON / (physicalconstants::K_B * electronTemperature())) - 1.);
            //if(retval < 0 || !isfinite(retval)) {
            //   retval = 0;
            //}
            //return retval;
         }
 
      };
 
      std::vector<Node> nodes;
 
      // Atmospheric height layers that are being integrated over
      constexpr static int numAtmosphereLevels = 20;
      struct AtmosphericLayer {
         Real altitude;
         Real nui;
         Real nue;
         Real density;
         Real depth; // integrated density from the top of the atmosphere
         Real pedersencoeff;
         Real hallcoeff;
         Real parallelcoeff;
      };
      std::array<AtmosphericLayer, numAtmosphereLevels> atmosphere;
 
      enum IonosphereSolverGaugeFixing { // Potential solver gauge fixing method
         None,     // No gauge fixing, solver won't converge well
         Pole,     // Fixing north pole (node 0) potential to zero
         Integral, // Fixing integral of potential to zero (unstable?)
         Equator   // Fixing all nodes within +-10 dgrees to zero
      } gaugeFixing;
 
      enum IonosphereIonizationModel { // Ionization production rate model
         Rees1963, // Rees (1963)
         Rees1989, // Rees (1989)
         SergienkoIvanov, // Sergienko & Ivanov (1993)
         Robinson2020, // Robinson et al (2020)
         Juusola2025,  // Juusola et al (2025)
         FixedSigma, // Simple, homogeneous, fixed SigmaP and SigmaH values, for testing
      } ionizationModel;
 
      // Ionisation production table
      std::array< std::array< std::array< Real, productionNumTemperatures >, productionNumAccEnergies >, numAtmosphereLevels > productionTable;
      Real lookupProductionValue(int heightindex, Real energy_keV, Real temperature_keV);
 
      MPI_Comm communicator = MPI_COMM_NULL; 
      int rank = -1;                         
      int writingRank;                       
      bool isCouplingInwards = true;         
      bool isCouplingOutwards = true;        
      FieldFunction dipoleField;             
      std::array<Real, 3> BGB;               
 
      std::map< std::array<Real, 3>, std::array<
         std::pair<int, Real>, 3> > vlasovGridCoupling; 
 
      void setDipoleField(const FieldFunction& dipole) {
         dipoleField = dipole;
      };
      void setConstantBackgroundField(const std::array<Real, 3> B) {
         BGB = B;
      }
      void readAtmosphericModelFile(const char* filename);
      void storeNodeB();
      void offset_FAC();                  
      void normalizeRadius(Node& n, Real R); 
      void updateConnectivity();          
      void updateIonosphereCommunicator(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid, fsgrids::technicalspan technical, FieldSolverGrid &fsgrid); 
      void initializeTetrahedron();       
      void initializeOctahedron();        
      void initializeIcosahedron();       
      void initializeSphericalFibonacci(int n); 
      void initializeGridFromFile(std::string path); 
      int32_t findElementNeighbour(uint32_t e, int n1, int n2);
      uint32_t findNodeAtCoordinates(std::array<Real,3> x); 
      void subdivideElement(uint32_t e);  
      void stitchRefinementInterfaces(); 
      void calculatePrecipitation(); 
      void calculateConductivityTensor(const Real F10_7, const Real recombAlpha, const Real backgroundIonisation, const bool refillTensorAtRestart=false); 
      Real interpolateUpmappedPotential(const std::array<Real, 3>& x); 
 
      // Conjugate Gradient solver functions
      void addMatrixDependency(uint node1, uint node2, Real coeff, bool transposed=false); 
      void addAllMatrixDependencies(uint nodeIndex);
      void initSolver(bool zeroOut=true);  
      iSolverReal Atimes(uint nodeIndex, int parameter, bool transpose=false); 
      Real Asolve(uint nodeIndex, int parameter, bool transpose=false); 
      void solve(
         int & iteration,
         int & nRestarts,
         Real & residual,
         Real & minPotentialN,
         Real & maxPotentialN,
         Real & minPotentialS,
         Real & maxPotentialS
      );
      void solveInternal(
         int & iteration,
         int & nRestarts,
         Real & residual,
         Real & minPotentialN,
         Real & maxPotentialN,
         Real & minPotentialS,
         Real & maxPotentialS
      );
 
      // Map field-aligned currents, density and temperature
      // down from the simulation boundary onto this grid
      void mapDownBoundaryData(fsgrids::perbspan perb,
                               fsgrids::constdperbspan dperb,
                               fsgrids::momentsspan moments,
                               fsgrids::technicalspan technical, FieldSolverGrid &fsgrid
      );
 
      // Returns the surface area of one element on the sphere
      Real elementArea(uint32_t elementIndex) {
         const std::array<Real, 3>& a = nodes[elements[elementIndex].corners[0]].x;
         const std::array<Real, 3>& b = nodes[elements[elementIndex].corners[1]].x;
         const std::array<Real, 3>& c = nodes[elements[elementIndex].corners[2]].x;
 
         // Two edges e1 = b-c,  e2 = c-a
         std::array<Real, 3> e1{b[0]-c[0], b[1]-c[1],b[2]-c[2]};
         std::array<Real, 3> e2{c[0]-a[0], c[1]-a[1],c[2]-a[2]};
         // Area vector A = cross(e1 e2)
         std::array<Real, 3> area{ e1[1]*e2[2] - e1[2]*e2[1],
                                   e1[2]*e2[0] - e1[0]*e2[2],
                                   e1[0]*e2[1] - e1[1]*e2[0]};
 
         return 0.5 * sqrt( area[0]*area[0] + area[1]*area[1] + area[2]*area[2] );
      }
 
      // Returns the projected surface area of one element, mapped up along the magnetic field to
      // the simulation boundary. If one of the nodes maps nowhere, returns 0.
      // Returns an oriented vector, which can be dotted with B
      std::array<Real, 3> mappedElementArea(uint32_t elementIndex) {
         const std::array<Real, 3>& a = nodes[elements[elementIndex].corners[0]].xMapped;
         const std::array<Real, 3>& b = nodes[elements[elementIndex].corners[1]].xMapped;
         const std::array<Real, 3>& c = nodes[elements[elementIndex].corners[2]].xMapped;
 
         // Check if any node maps to zero
         if( sqrt( a[0]*a[0] + a[1]*a[1] + a[2]*a[2] ) == 0 ||
               sqrt( b[0]*b[0] + b[1]*b[1] + b[2]*b[2] ) == 0 ||
               sqrt( c[0]*c[0] + c[1]*c[1] + c[2]*c[2] ) == 0) {
 
            return {0,0,0};
         }
 
         // Two edges e1 = b-c,  e2 = c-a
         std::array<Real, 3> e1{b[0]-c[0], b[1]-c[1],b[2]-c[2]};
         std::array<Real, 3> e2{c[0]-a[0], c[1]-a[1],c[2]-a[2]};
         // Area vector A = cross(e1 e2)
         const Real HALF = 0.5;
         const Real THIRD = 1./3.;
         std::array<Real, 3> area{ HALF * (e1[1]*e2[2] - e1[2]*e2[1]),
                                   HALF * (e1[2]*e2[0] - e1[0]*e2[2]),
                                   HALF * (e1[0]*e2[1] - e1[1]*e2[0])};
 
         // By definition, the area is oriented outwards, so if dot(r,A) < 0, flip it.
         std::array<Real, 3> r{
            (a[0]+b[0]+c[0]) * THIRD,
            (a[1]+b[1]+c[1]) * THIRD,
            (a[2]+b[2]+c[2]) * THIRD};
         if(area[0]*r[0] + area[1]*r[1] + area[2] *r[2] < 0) {
            area[0]*=-1.;
            area[1]*=-1.;
            area[2]*=-1.;
         }
         return area;
      }
 
      Eigen::Vector3d elementBarycentre(uint32_t el) {
         Eigen::Vector3d barycentre(0,0,0);
 
         Element& element = elements[el];
         for(uint i=0; i<3; i++) {
            Eigen::Vector3d corner(nodes[element.corners[i]].x.data());
 
            barycentre += corner;
         }
         barycentre /= 3.;
 
         return barycentre;
      }
 
      Eigen::Vector3d elementCircumcentre(uint el) {
         Eigen::Vector3d circumcentre(0,0,0);
 
         SphericalTriGrid::Element& element = elements[el];
         uint corner1 = element.corners[0];
         uint corner2 = element.corners[1];
         uint corner3 = element.corners[2];
 
         Eigen::Vector3d a(nodes[corner1].x.data());
         Eigen::Vector3d b(nodes[corner2].x.data());
         Eigen::Vector3d c(nodes[corner3].x.data());
 
         Eigen::Vector3d edge1 = b - a;
         Eigen::Vector3d edge2 = c - a;
 
         Eigen::Vector3d edge1Mid = a + edge1 / 2.;
         Eigen::Vector3d edge2Mid = a + edge2 / 2.;
 
         Eigen::Vector3d normal = edge1.cross(edge2).normalized();
 
         if(normal.dot(a) < 0) {
            normal *= -1.;
         }
 
         Eigen::Vector3d edge1Perpendicular = normal.cross(edge1).normalized();
         Eigen::Vector3d edge2Perpendicular = normal.cross(edge2).normalized();
 
         Eigen::Matrix<Real, 3, 2> A;
         A.col(0) = edge1Perpendicular;
         A.col(1) = - edge2Perpendicular;
         Eigen::Vector3d bVec = edge2Mid - edge1Mid;
         Eigen::Vector2d t = A.colPivHouseholderQr().solve(bVec);
         // Verify that the solution is correct
         //Eigen::Vector3d residual = A * t - bVec;
         //if (residual.norm() > 1e-6) {
         //   cerr << "Circumcentre calculation failed, residual: " << residual.norm() << endl;
         //}
 
         Eigen::Vector3d intersection = edge1Mid + t(0) * edge1Perpendicular;
         Eigen::Vector3d intersection2 = edge2Mid + t(1) * edge2Perpendicular;
         if((intersection - intersection2).norm() > 1e-6) {
            cerr << "Circumcentre calculation failed, intersection points do not match: "
                 << (intersection - intersection2).norm() << endl;
         }
         circumcentre = intersection;
 
         return circumcentre;
      }
 
      Eigen::Vector3d elementNormal(uint32_t el) {
         Eigen::Vector3d normal(0,0,0);
 
         SphericalTriGrid::Element& element = elements[el];
         uint32_t corner1 = element.corners[0];
         uint32_t corner2 = element.corners[1];
         uint32_t corner3 = element.corners[2];
 
         Eigen::Vector3d a(nodes[corner1].x.data());
         Eigen::Vector3d b(nodes[corner2].x.data());
         Eigen::Vector3d c(nodes[corner3].x.data());
 
         Eigen::Vector3d edge1 = b - a;
         Eigen::Vector3d edge2 = c - a;
 
         normal = edge1.cross(edge2);
 
         normal.normalized();
 
         if(normal.dot(elementCircumcentre(el)) < 0) {
            normal *= -1.;
         }
 
         return normal;
      }
 
      Eigen::Vector3d commonEdgeMidpoint(uint32_t el1, uint32_t el2) {
         SphericalTriGrid::Element& element1 = elements[el1];
         SphericalTriGrid::Element& element2 = elements[el2];
 
         // Get common edge to these two elements
         for(uint i=0; i<3; i++) {
            if(element1.corners[i] == element2.corners[0] ||
               element1.corners[i] == element2.corners[1] ||
               element1.corners[i] == element2.corners[2]) {
               for(uint j=0; j<3; j++) {
                  if(i != j && (element1.corners[j] == element2.corners[0] ||
                              element1.corners[j] == element2.corners[1] ||
                              element1.corners[j] == element2.corners[2])) {
 
                        uint corner1 = element1.corners[i];
                        uint corner2 = element1.corners[j];
 
                        Eigen::Vector3d a(nodes[corner1].x.data());
                        Eigen::Vector3d b(nodes[corner2].x.data());
 
                        Eigen::Vector3d midpoint = (a + b) / 2.;
                        return midpoint;
                  }
               }
            }
         }
 
        // We should never end up here.
        return {0,0,0};
      }
 
      Real dualPolygonArea(uint gridNode){
         Real A = 0.;
 
         for(uint i = 0; i < nodes[gridNode].numTouchingElements; i++){
            uint32_t gridEl = nodes[gridNode].touchingElements[i];
            Eigen::Vector3d nodePosition(nodes[gridNode].x.data());
 
            SphericalTriGrid::Element& element = elements[gridEl];
 
            int localC=0,localI=0,localJ=0;
            for(int c=0; c < 3; c++) {
               if(element.corners[c] == gridNode) {
                  localC = c;
                  localI = (c+1)%3;
                  localJ = (c+2)%3;
                  break;
               }
            }
 
            uint otherElementi = findElementNeighbour(gridEl, localC, localI);
            uint otherElementj = findElementNeighbour(gridEl, localC, localJ);
 
            Eigen::Vector3d midpointi = commonEdgeMidpoint(gridEl, otherElementi);
            Eigen::Vector3d midpointj = commonEdgeMidpoint(gridEl, otherElementj);
 
            Eigen::Vector3d circumcentre = elementCircumcentre(gridEl);
 
            Real heighti = (nodePosition - midpointi).norm();
            Real heightj = (nodePosition - midpointj).norm();
 
            Real basei = (circumcentre - midpointi).norm();
            Real basej = (circumcentre - midpointj).norm();
 
            A += (0.5 * basei * heighti) + (0.5 * basej * heightj);
         }
 
         return A;
      }
 
      Real areaInDualPolygon(uint gridNode, uint gridElem) {
         Real A = 0.;
 
         Eigen::Vector3d nodePosition(nodes[gridNode].x.data());
         SphericalTriGrid::Element& element = elements[gridElem];
 
         int localC=0,localI=0,localJ=0;
         for(int c=0; c < 3; c++) {
            if(element.corners[c] == gridNode) {
               localC = c;
               localI = (c+1)%3;
               localJ = (c+2)%3;
               break;
            }
         }
 
         uint otherElementi = findElementNeighbour(gridElem, localC, localI);
         uint otherElementj = findElementNeighbour(gridElem, localC, localJ);
 
         Eigen::Vector3d midpointi = commonEdgeMidpoint(gridElem, otherElementi);
         Eigen::Vector3d midpointj = commonEdgeMidpoint(gridElem, otherElementj);
 
         Eigen::Vector3d circumcentre = elementCircumcentre(gridElem);
 
         Real heighti = (nodePosition - midpointi).norm();
         Real heightj = (nodePosition - midpointj).norm();
 
         Real basei = (circumcentre - midpointi).norm();
         Real basej = (circumcentre - midpointj).norm();
 
         A += (0.5 * basei * heighti) + (0.5 * basej * heightj);
 
         return A;
      }
 
      Real nodeNeighbourArea(uint32_t nodeIndex) { // Summed area of all touching elements
 
         Node& n = nodes[nodeIndex];
         Real area=0;
 
         for(uint i=0; i<n.numTouchingElements; i++) {
            area += elementArea(n.touchingElements[i]);
         }
         return area;
      }
 
      // Calculate neighbor's Barycentre and dual polygon - edge - intersection point.
      std::tuple<Eigen::Vector3d, Eigen::Vector3d> connectingSegmentLengths(uint32_t el1, uint32_t el2) {
         SphericalTriGrid::Element& element1 = elements[el1];
         SphericalTriGrid::Element& element2 = elements[el2];
 
         Eigen::Vector3d barycentre1 = elementBarycentre(el1);
         Eigen::Vector3d barycentre2 = elementBarycentre(el2);
 
         // Get common edge to these two elements
         for(uint i=0; i<3; i++) {
            if(element1.corners[i] == element2.corners[0] ||
                  element1.corners[i] == element2.corners[1] ||
                  element1.corners[i] == element2.corners[2]) {
               for(uint j=0; j<3; j++) {
                  if(i != j && (element1.corners[j] == element2.corners[0] ||
                           element1.corners[j] == element2.corners[1] ||
                           element1.corners[j] == element2.corners[2])) {
 
                     uint corner1 = element1.corners[i];
                     uint corner2 = element1.corners[j];
 
                     Eigen::Vector3d normal1 = elementNormal(el1);
                     Eigen::Vector3d normal2 = elementNormal(el2);
 
                     Eigen::Vector3d rotatedBarycentre2 =  Eigen::Vector3d(nodes[corner1].x.data()) +
                        Eigen::Quaternion<Real>::FromTwoVectors(normal2, normal1).toRotationMatrix() *
                        (barycentre2 - Eigen::Vector3d(nodes[corner1].x.data()));
 
                     Eigen::Vector3d corner1Position(nodes[corner1].x.data());
                     Eigen::Vector3d corner2Position(nodes[corner2].x.data());
 
                     Eigen::Vector3d barycentre1ToBarycentre2 = (rotatedBarycentre2 - barycentre1).normalized();
                     Eigen::Vector3d corner1ToCorner2 = (corner2Position - corner1Position).normalized();
 
                     // Get intersection of line between barycenters and line between corners
                     Eigen::Matrix<double, 3, 2> A;
                     A.col(0) = barycentre1ToBarycentre2;
                     A.col(1) = - corner1ToCorner2;
                     Eigen::Vector3d b = corner1Position - barycentre1;
                     Eigen::Vector2d t = A.colPivHouseholderQr().solve(b);
                     Eigen::Vector3d intersection = barycentre1 + t(0) * barycentre1ToBarycentre2;
 
                     return std::make_tuple(barycentre2, intersection);
                  }
               }
            }
         }
 
         // Not found, something went bananas.
         abort();
      }
 
      std::array<Real,3> computeGradT(const std::array<Real, 3>& a, const std::array<Real, 3>& b, const std::array<Real, 3>& c);
      std::array<Real, 9> sigmaAverage(uint elementIndex);
      double elementIntegral(uint elementIndex, int i, int j, bool transpose = false);
 
   };
 
   extern SphericalTriGrid ionosphereGrid;

   class Ionosphere: public SysBoundaryCondition {
   public:
      Ionosphere();
      virtual ~Ionosphere();
 
      static void addParameters();
      virtual void getParameters() override;
 
      virtual void initSysBoundary(creal& t, Project& project) override;
      virtual void assignSysBoundary(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                     fsgrids::technicalspan technical, FieldSolverGrid& fsgrid) override;
      virtual void applyInitialState(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                     fsgrids::technicalspan technical, FieldSolverGrid& fsgrid,
                                     fsgrids::perbspan perb,
                                     fsgrids::bgbspan bgb,
                                     Project& project) override;
      virtual Real fieldSolverBoundaryCondMagneticField(fsgrids::perbspan b,
                                                        fsgrids::constbgbspan bgb,
                                                        fsgrids::consttechnicalspan technical,
                                                        const std::array<Real, 3>& gridSpacing,
                                                        const std::array<fsgrid::FsSize_t, 3>& globalCoordinates,
                                                        const fsgrid::FsStencil& stencil, cuint component) override;
      virtual void fieldSolverBoundaryCondElectricField(fsgrids::efieldspan e,
                                                        const fsgrid::FsStencil& stencil, cuint component) override;
      virtual void fieldSolverBoundaryCondHallElectricField(fsgrids::ehallspan ehall,
                                                            const fsgrid::FsStencil& stencil, cuint component) override;
      virtual void fieldSolverBoundaryCondGradPeElectricField(fsgrids::egradpespan EGradPe,
                                                              const fsgrid::FsStencil& stencil, cuint component) override;
      virtual void fieldSolverBoundaryCondDerivatives(fsgrids::dperbspan dperb,
                                                      fsgrids::dmomentsspan dmoments,
                                                      const fsgrid::FsStencil& stencil, cuint RKCase, cuint component) override;
      virtual void fieldSolverBoundaryCondBVOLDerivatives(fsgrids::volspan vols,
                                                          const fsgrid::FsStencil& stencil, cuint component) override;
      // Compute and store the EXB drift into the cell's BULKV_FORCING_X/Y/Z fields
      virtual void mapCellPotentialAndGetEXBDrift(std::array<Real, CellParams::N_SPATIAL_CELL_PARAMS>& cellParams) override;
      virtual void vlasovBoundaryCondition(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                           const CellID& cellID, const uint popID,
                                           const bool calculate_V_moments) override;
      virtual void updateState(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                               fsgrids::technicalspan technical, FieldSolverGrid& fsgrid,
                               fsgrids::perbspan perb,
                               fsgrids::bgbspan bgb, creal t) override;
 
      virtual void getFaces(bool* faces) override;
      virtual std::string getName() const override;
      virtual uint getIndex() const override;
      static Real radius; 
      static std::vector<IonosphereSpeciesParameters> speciesParams;
 
      // Parameters of the ionosphere model
      static Real innerRadius; 
      static bool useEigenSolver;
      static int solverMaxIterations; 
      static Real solverRelativeL2ConvergenceThreshold; 
      static int solverMaxFailureCount;
      static Real solverMaxErrorGrowthFactor;
      static bool solverPreconditioning; 
      static bool solverUseMinimumResidualVariant; 
      static bool solverToggleMinimumResidualVariant; 
      static Real shieldingLatitude; 
      static Real ridleyParallelConductivity; 
 
      // TODO: Make these parameters of the IonosphereGrid
      static Real recombAlpha; 
      static Real F10_7; 
      static Real backgroundIonisation; 
      static Real downmapRadius; 
      static Real downmapSamplingWidth; 
      static enum downmapSamplingMode { // How to sample possibly under-resolved FACs at the downmap radius
         Pointwise,  // Just sample the FAC at the downmapping point
         Boxcar27   // 27-point boxcar, samples a cube of +-downmapSamplingWidth*fsgrid.dx/2 from the downmapping point
      } downmapFACsamplingMode;
      
      static Real unmappedNodeRho; 
      static Real unmappedNodeTe; 
      static Real couplingTimescale; 
      static Real couplingInterval; 
      static int solveCount; 
      static enum IonosphereConductivityModel { // How should the conductivity tensor be assembled?
         GUMICS,   // Like GUMICS-5 does it? (Only SigmaH and SigmaP, B perp to surface)
         Ridley,   // Or like the Ridley 2004 paper (with 1000 mho longitudinal conductivity)
         Koskinen  // Like Koskinen's 2001 "Physics of Space Storms" book suggests
      } conductivityModel;
      static Real fixedSigmaP; 
      static Real fixedSigmaH;
 
      void generateTemplateCell(Project &project);
      void setCellFromTemplate(SpatialCell* cell,const uint popID);
 
      std::array<Real, 3> fieldSolverGetNormalDirection(
         fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
         cint i,
         cint j,
         cint k
      );
 
      Real center[3]; 
      uint geometry; 
 
 
      std::string baseShape; 
      std::string path; 
      int fibonacciNodeNum;  
      Real earthAngularVelocity; 
      Real plasmapauseL; 
      std::string atmosphericModelFile; 
 
      // Boundaries of refinement latitude bands
      std::vector<Real> refineMinLatitudes;
      std::vector<Real> refineMaxLatitudes;
 
      spatial_cell::SpatialCell templateCell;
   };
}
 
#endif