Started Electril Sail simulation project

- fixed bug in MultiPeak adaptive sampling - added 2D solvers to Poisson
philowzx · Mar 5, 2015 · 55efa37 · 55efa37
1 parent ec66e56
commit 55efa37
Show file tree

Hide file tree

Showing 16 changed files with 786 additions and 53 deletions.
diff --git a/Makefile b/Makefile
@@ -138,6 +138,7 @@ DEPS_PROJECTS =	projects/project.h projects/project.cpp \
 		projects/Diffusion/Diffusion.h projects/Diffusion/Diffusion.cpp \
 		projects/Dispersion/Dispersion.h projects/Dispersion/Dispersion.cpp \
 		projects/Distributions/Distributions.h projects/Distributions/Distributions.cpp \
+		projects/ElectricSail/electric_sail.h projects/ElectricSail/electric_sail.cpp \
 		projects/Firehose/Firehose.h projects/Firehose/Firehose.cpp \
 		projects/Flowthrough/Flowthrough.h projects/Flowthrough/Flowthrough.cpp \
 		projects/Fluctuations/Fluctuations.h projects/Fluctuations/Fluctuations.cpp \
@@ -189,7 +190,7 @@ OBJS = 	version.o memoryallocation.o backgroundfield.o quadr.o dipole.o linedipo
 	donotcompute.o ionosphere.o outflow.o setbyuser.o setmaxwellian.o \
 	sysboundary.o sysboundarycondition.o \
 	project.o projectTriAxisSearch.o \
-	Alfven.o Diffusion.o Dispersion.o Distributions.o Firehose.o Flowthrough.o Fluctuations.o Harris.o KHB.o Larmor.o \
+	Alfven.o Diffusion.o Dispersion.o Distributions.o electric_sail.o Firehose.o Flowthrough.o Fluctuations.o Harris.o KHB.o Larmor.o \
 	Magnetosphere.o MultiPeak.o VelocityBox.o Riemann1.o Shock.o Template.o test_fp.o testHall.o test_trans.o \
 	verificationLarmor.o Shocktest.o grid.o ioread.o iowrite.o vlasiator.o logger.o muxml.o \
 	common.o parameters.o readparameters.o spatial_cell.o \
@@ -296,6 +297,9 @@ Dispersion.o: ${DEPS_COMMON} projects/Dispersion/Dispersion.h projects/Dispersio
 Distributions.o: ${DEPS_COMMON} projects/Distributions/Distributions.h projects/Distributions/Distributions.cpp
 	${CMP} ${CXXFLAGS} ${FLAGS} ${MATHFLAGS} -c projects/Distributions/Distributions.cpp ${INC_DCCRG} ${INC_ZOLTAN} ${INC_BOOST} ${INC_EIGEN}
 
+electric_sail.o: ${DEPS_COMMON} projects/ElectricSail/electric_sail.h projects/ElectricSail/electric_sail.cpp
+	${CMP} ${CXXFLAGS} ${FLAGS} ${MATHFLAGS} -c projects/ElectricSail/electric_sail.cpp ${INC_DCCRG} ${INC_ZOLTAN} ${INC_BOOST} ${INC_EIGEN}
+
 Firehose.o: ${DEPS_COMMON} projects/Firehose/Firehose.h projects/Firehose/Firehose.cpp
 	${CMP} ${CXXFLAGS} ${FLAGS} ${MATHFLAGS} -c projects/Firehose/Firehose.cpp ${INC_DCCRG} ${INC_ZOLTAN} ${INC_BOOST} ${INC_EIGEN}
 

diff --git a/common.h b/common.h
@@ -220,9 +220,12 @@ namespace CellParams {
       LBWEIGHTCOUNTER,    /*!< Counter for storing compute time weights needed by the load balancing**/
       ACCSUBCYCLES,        /*!< number of subcyles for each cell*/
       ISCELLSAVINGF,      /*!< Value telling whether a cell is saving its distribution function when partial f data is written out. */
-      PHI,      /*!< Electrostatic potential.*/
-      PHI_TMP,  /*!< Temporary electrostatic potential.*/
-      RHOQ_TOT, /*!< Total charge density, summed over all particle populations.*/
+      PHI,        /*!< Electrostatic potential.*/
+      PHI_TMP,    /*!< Temporary electrostatic potential.*/
+      RHOQ_TOT,   /*!< Total charge density, summed over all particle populations.*/
+      BGEXVOL,    /*!< Background electric field averaged over spatial cell, x-component.*/
+      BGEYVOL,    /*!< Background electric field averaged over spatial cell, y-component.*/
+      BGEZVOL,    /*!< Background electric field averaged over spatial cell, z-component.*/
       N_SPATIAL_CELL_PARAMS
    };
 }

diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp
@@ -1,18 +1,7 @@
 /*
 This file is part of Vlasiator.
 
-Copyright 2010, 2011, 2012, 2013 Finnish Meteorological Institute
-
-
-
-
-
-
-
-
-
-
-
+Copyright 2010-2013,2015 Finnish Meteorological Institute
 
 */
 
@@ -222,6 +211,9 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
       if (*it == "Potential") {
          outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("poisson/potential",CellParams::PHI,1));
       }
+      if (*it == "BackgroundVolE") {
+         outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("poisson/BGE_vol",CellParams::BGEXVOL,3));
+      }
       if (*it == "ChargeDensity") {
          outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("poisson/rho_q",CellParams::RHOQ_TOT,1));
       }

diff --git a/poisson_solver/poisson_solver.cpp b/poisson_solver/poisson_solver.cpp
@@ -1,7 +1,7 @@
 /* This file is part of Vlasiator.
  * Copyright 2015 Finnish Meteorological Institute.
  * 
- * File:   poisson_solver.h
+ * File:   poisson_solver.cpp
  * Author: sandroos
  *
  * Created on January 14, 2015, 1:42 PM
@@ -33,6 +33,7 @@ namespace poisson {
    int Poisson::PHI = CellParams::PHI;
    ObjectFactory<PoissonSolver> Poisson::solvers;
    PoissonSolver* Poisson::solver = NULL;
+   bool Poisson::is2D = false;
    string Poisson::solverName;
    uint Poisson::maxIterations;
    Real Poisson::minRelativePotentialChange;
@@ -63,6 +64,99 @@ namespace poisson {
 
    bool PoissonSolver::finalize() {return true;}
 
+   /** Calculate total charge density on given spatial cells.
+    * @param mpiGrid Parallel grid library.
+    * @param cells List of spatial cells.
+    * @return If true, charge densities were successfully calculated.*/
+   bool PoissonSolver::calculateChargeDensity(spatial_cell::SpatialCell* cell) {
+      phiprof::start("Charge Density");
+      bool success = true;
+
+      Real rho_q = 0.0;
+      #pragma omp parallel reduction (+:rho_q)
+      {
+         // Iterate all particle species
+         for (int popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
+            vmesh::VelocityBlockContainer<vmesh::LocalID>& blockContainer = cell->get_velocity_blocks(popID);
+            if (blockContainer.size() == 0) continue;
+
+            const Real charge       = getObjectWrapper().particleSpecies[popID].charge;
+            const Realf* data       = blockContainer.getData();
+            const Real* blockParams = blockContainer.getParameters();
+
+            // Sum charge density over all phase-space cells
+            #pragma omp for
+            for (vmesh::LocalID blockLID=0; blockLID<blockContainer.size(); ++blockLID) {
+               Real sum = 0.0;
+               for (int i=0; i<WID3; ++i) sum += data[blockLID*WID3+i];
+
+               const Real DV3 
+                  = blockParams[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS+BlockParams::DVX]
+                  * blockParams[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS+BlockParams::DVY]
+                  * blockParams[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS+BlockParams::DVZ];
+               rho_q += charge*sum/DV3;
+            }
+         }
+      }
+      cell->parameters[CellParams::RHOQ_TOT] = rho_q;
+
+      phiprof::stop("Charge Density");
+      return success;
+   }
+
+   bool PoissonSolver::calculateElectrostaticField2D(const std::vector<poisson::CellCache3D>& cells) {
+      phiprof::start("Electrostatic E");
+
+      #pragma omp parallel for
+      for (size_t c=0; c<cells.size(); ++c) {
+         const Real DX = 2*cells[c].parameters[0][CellParams::DX];
+         const Real DY = 2*cells[c].parameters[0][CellParams::DY];
+
+         const Real phi_01 = cells[c].parameters[1][CellParams::PHI]; // -x neighbor
+         const Real phi_21 = cells[c].parameters[2][CellParams::PHI]; // +x neighbor
+         const Real phi_10 = cells[c].parameters[3][CellParams::PHI]; // -y neighbor
+         const Real phi_12 = cells[c].parameters[4][CellParams::PHI]; // +y neighbor
+
+         cells[c].parameters[0][CellParams::EXVOL] -= (phi_21-phi_01)/DX;
+         cells[c].parameters[0][CellParams::EYVOL] -= (phi_12-phi_10)/DY;
+      }
+
+      phiprof::stop("Electrostatic E");
+      return true;
+   }
+
+   bool PoissonSolver::calculateElectrostaticField3D(const std::vector<poisson::CellCache3D>& cells) {
+      phiprof::start("Electrostatic E");
+
+      #pragma omp parallel for
+      for (size_t c=0; c<cells.size(); ++c) {
+         const Real DX = 2*cells[c].parameters[0][CellParams::DX];
+         const Real DY = 2*cells[c].parameters[0][CellParams::DY];
+         const Real DZ = 2*cells[c].parameters[0][CellParams::DZ];
+
+         const Real phi_011 = cells[c].parameters[1][CellParams::PHI]; // -x neighbor
+         const Real phi_211 = cells[c].parameters[2][CellParams::PHI]; // +x neighbor
+         const Real phi_101 = cells[c].parameters[3][CellParams::PHI]; // -y neighbor
+         const Real phi_121 = cells[c].parameters[4][CellParams::PHI]; // +y neighbor
+         const Real phi_110 = cells[c].parameters[5][CellParams::PHI]; // -z neighbor
+         const Real phi_112 = cells[c].parameters[6][CellParams::PHI]; // +z neighbor
+
+         cells[c].parameters[0][CellParams::EXVOL] -= (phi_211-phi_011)/DX;
+         cells[c].parameters[0][CellParams::EYVOL] -= (phi_121-phi_101)/DY;
+         cells[c].parameters[0][CellParams::EZVOL] -= (phi_112-phi_110)/DZ;
+      }
+
+      phiprof::stop("Electrostatic E");
+      return true;
+   }
+
+   /** Estimate the error in the numerical solution of the electrostatic potential.
+    * The error is calculated as the difference of new and old potential, divided 
+    * by the new potential value. This function works for both 2D and 3D solver.
+    * This function must be called simultaneously by all MPI processes.
+    * @param mpiGrid Parallel grid.
+    * @return The relative error in Poisson equation solution. The return value is 
+    * the same at all MPI processes.*/
    Real PoissonSolver::error(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) {
       phiprof::start("Potential Change");
 
@@ -215,10 +309,16 @@ namespace poisson {
       }
 
       // Set up the initial state unless the simulation was restarted
-      vector<CellID> local_cells = mpiGrid.get_cells();
+      if (Parameters::isRestart == true) return success;
+
+      const vector<CellID>& local_cells = getLocalCells();
       for (size_t c=0; c<local_cells.size(); ++c) {
          spatial_cell::SpatialCell* cell = mpiGrid[local_cells[c]];
          cell->parameters[CellParams::PHI] = 0;
+         cell->parameters[CellParams::PHI_TMP] = 0;
+         cell->parameters[CellParams::EXVOL] = cell->parameters[CellParams::BGEXVOL];
+         cell->parameters[CellParams::EYVOL] = cell->parameters[CellParams::BGEYVOL];
+         cell->parameters[CellParams::EZVOL] = cell->parameters[CellParams::BGEZVOL];
       }
 
       return success;
@@ -236,7 +336,8 @@ namespace poisson {
 
    bool solve(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) {
       phiprof::start("Poisson Solver (Total)");
-
+      bool success = true;
+
       // If mesh partitioning has changed, recalculate spatial 
       // cell parameters pointer cache:
       if (Parameters::meshRepartitioned == true) {
@@ -245,10 +346,15 @@ namespace poisson {
 	 phiprof::stop("Cache Cell Parameters");
       }
 
-      bool success = true;
-      if (Poisson::solver != NULL) {
+      // Solve Poisson equation
+      if (success == true) if (Poisson::solver != NULL) {
          if (Poisson::solver->solve(mpiGrid) == false) success = false;
       }
+
+      // Add electrostatic electric field to volume-averaged E
+      if (success == true) if (Poisson::solver != NULL) {
+         if (Poisson::solver->calculateElectrostaticField(mpiGrid) == false) success = false;
+      }
 
       phiprof::stop("Poisson Solver (Total)");
       return success;

diff --git a/poisson_solver/poisson_solver.h b/poisson_solver/poisson_solver.h
@@ -19,11 +19,13 @@
 namespace poisson {
 
     struct CellCache2D {
+       spatial_cell::SpatialCell* cell;
        Real* parameters[5];
        Real*& operator[](const int& i) {return parameters[i];}
     };
 
     struct CellCache3D {
+       spatial_cell::SpatialCell* cell;
        Real* parameters[7];
        Real*& operator[](const int& i) {return parameters[i];}
     };
@@ -33,10 +35,26 @@ namespace poisson {
        PoissonSolver();
        virtual ~PoissonSolver();
 
+       // ***** DECLARATIONS OF VIRTUAL MEMBER FUNCTIONS ***** //
+
+       virtual bool calculateChargeDensity(spatial_cell::SpatialCell* cell);
+       virtual bool calculateElectrostaticField2D(const std::vector<poisson::CellCache3D>& cells);
+       virtual bool calculateElectrostaticField3D(const std::vector<poisson::CellCache3D>& cells);
        virtual Real error(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
        virtual Real error3D(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
        virtual bool initialize();
        virtual bool finalize();
+
+       // ***** DECLARATIONS OF PURE VIRTUAL MEMBER FUNCTIONS ***** //
+
+       /** Calculate electric field on all non-system boundary spatial cells.
+        * @param mpiGrid Parallel grid library.
+        * @return If true, electric field was calculated successfully.*/
+       virtual bool calculateElectrostaticField(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) = 0;
+
+       /** Solve Poisson equation on all non-system boundary spatial cells.
+        * @param mpiGrid Parallel grid library.
+        * @return If true, Poisson equation was successfully solved.*/
        virtual bool solve(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) = 0;
 
     protected:
@@ -52,6 +70,9 @@ namespace poisson {
       static PoissonSolver* solver;                /**< Poisson solver used in the simulation.*/
       static std::string solverName;               /**< Name of the Poisson solver in use.*/
 
+      static bool is2D;                            /**< If true, then system is two-dimensional, i.e., 
+                                                    * electrostatic potential and electric field is solved 
+                                                    * in xy-plane.*/
       static uint maxIterations;                   /**< Maximum number of iterations allowed, only 
                                                     * has an effect on iterative solvers.*/
       static Real minRelativePotentialChange;      /**< Iterative solvers keep on iterating the solution 

diff --git a/poisson_solver/poisson_solver_jacobi.cpp b/poisson_solver/poisson_solver_jacobi.cpp
@@ -39,6 +39,11 @@ namespace poisson {
       return success;
    }
 
+   bool PoissonSolverJacobi::calculateElectrostaticField(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid) {
+      #warning Jacobi solver does not calculate electric field
+      return false;
+   }
+
    void PoissonSolverJacobi::evaluate(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
                                    const std::vector<CellID>& cells) {      
       for (size_t c=0; c<cells.size(); ++c) {

diff --git a/poisson_solver/poisson_solver_jacobi.h b/poisson_solver/poisson_solver_jacobi.h
@@ -21,6 +21,7 @@ namespace poisson {
         PoissonSolverJacobi();
         ~PoissonSolverJacobi();
 
+        bool calculateElectrostaticField(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);
         bool initialize();
         bool finalize();
         bool solve(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid);