sysboundarycondition.h

Go to the documentation of this file.

/*
 * This file is part of Vlasiator.
 * Copyright 2010-2016 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 SYSBOUNDARYCONDITION_H
#define SYSBOUNDARYCONDITION_H
 
#include <dccrg.hpp>
#include <dccrg_cartesian_geometry.hpp>
#include <fsgrid.hpp>
 
#include <vector>
#include "../common.h"
#include "../definitions.h"
#include "../spatial_cells/spatial_cell_wrapper.hpp"
#include "../projects/project.h"
 
using namespace spatial_cell;
using namespace projects;
 
namespace SBC {
   class SysBoundaryCondition {
      public:
         SysBoundaryCondition();
         virtual ~SysBoundaryCondition();
 
         static void addParameters();
         virtual void getParameters() {
            std::cerr << "ERROR: base class SysBoundaryCondition::getParameters called!" << std::endl;
         }
         virtual void generateTemplateCell() {
            std::cerr << "ERROR: base class SysBoundaryCondition::generateTemplateCells called!" << std::endl;
         }
 
         virtual void initSysBoundary(
            creal& t,
            Project &project
         )=0;
         virtual void assignSysBoundary(dccrg::Dccrg<SpatialCell,
                                        dccrg::Cartesian_Geometry>& mpiGrid,
                                        fsgrids::technicalspan technical, FieldSolverGrid &fsgrid)=0;
         virtual void applyInitialState(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                        fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
                                        fsgrids::perbspan perb,
                                        fsgrids::bgbspan bgb,
                                        Project& project) = 0;
         virtual void updateState(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                  fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
                                  fsgrids::perbspan perb,
                                  fsgrids::bgbspan bgb, creal t) = 0;
         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) = 0;
         virtual void fieldSolverBoundaryCondElectricField(fsgrids::efieldspan e,
                                                           const fsgrid::FsStencil& stencil, cuint component) = 0;
         virtual void
         fieldSolverBoundaryCondHallElectricField(fsgrids::ehallspan ehall,
                                                  const fsgrid::FsStencil& stencil, cuint component) = 0;
         virtual void
         fieldSolverBoundaryCondGradPeElectricField(fsgrids::egradpespan EGradPe,
                                                    const fsgrid::FsStencil& stencil, cuint component) = 0;
         virtual void
         fieldSolverBoundaryCondDerivatives(fsgrids::dperbspan dperb,
                                            fsgrids::dmomentsspan dmoments,
                                            const fsgrid::FsStencil& stencil, cuint RKCase, cuint component) = 0;
         virtual void
         fieldSolverBoundaryCondBVOLDerivatives(fsgrids::volspan vols,
                                                const fsgrid::FsStencil& stencil, cuint component) = 0;
         static void setCellDerivativesToZero(fsgrids::dperbspan dperb,
                                              fsgrids::dmomentsspan dmoments,
                                              const fsgrid::FsStencil& stencil, cuint component);
         static void setCellBVOLDerivativesToZero(fsgrids::volspan vols,
                                                  const fsgrid::FsStencil& stencil, cuint component);
 
         virtual void mapCellPotentialAndGetEXBDrift(
            std::array<Real, CellParams::N_SPATIAL_CELL_PARAMS>& cellParams
         );

        virtual void vlasovBoundaryCondition(
            dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            const CellID& cellID,
            const uint popID,
            const bool calculate_V_moments
        )=0;
 
        virtual void setupL2OutflowAtRestart(
            dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid
        ) {
            std::cerr << "ERROR: base class SysBoundaryCondition::setupL2OutflowAtRestart called!" << std::endl;
        }
 
 

         virtual void getFaces(bool *faces) = 0;
         virtual std::string getName() const {
            std::cerr << "ERROR: base class SysBoundaryCondition::getName called!" << std::endl;
            return "ERROR";
         }
         virtual uint getIndex() const = 0; // {
//            std::cerr << "ERROR: base class SysBoundaryCondition::getIndex called!" << std::endl;
//            return sysboundarytype::N_SYSBOUNDARY_CONDITIONS;
//         }
         uint getPrecedence() const;
         bool isDynamic() const;
 
         bool updateSysBoundaryConditionsAfterLoadBalance(
            const dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            const std::vector<CellID> & local_cells_on_boundary
         );
         bool doApplyUponRestart() const;
         void setPeriodicity(
            std::array<bool, 3> isFacePeriodic
         );
      protected:
         void determineFace(
            bool* isThisCellOnAFace,
            const creal x,const  creal y,const creal z,
            const creal dx,const creal dy,const creal dz,
            const bool excludeSlicesAndPeriodicDimensions = false
         ) const;
         void determineFace(
            std::array<bool, 6> &isThisCellOnAFace,
            const dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            CellID id,
            const bool excludeSlicesAndPeriodicDimensions = false
         );
         void copyCellData(
            const SpatialCell *from,
            SpatialCell *to,
            const bool copyMomentsOnly,
            const uint popID,
            const bool copy_V_moments
         );
         std::array<SpatialCell*,27> & getFlowtoCells(
               const CellID& cellID
         );
 
         std::array<Realf*,27> getFlowtoCellsBlock(
               const std::array<SpatialCell*,27> flowtoCells,
               const vmesh::GlobalID blockGID,
               const uint popID
         );
 

      inline int nbrID(const int i, const int j, const int k){
         return (k+1)*9 + (j+1)*3 + i + 1;
      }
 
         void vlasovBoundaryCopyFromTheClosestNbr(
            dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            const CellID& cellID,
            const bool& copyMomentsOnly,
            const uint popID,
            const bool calculate_V_moments
         );
         void vlasovBoundaryCopyFromTheClosestL1OutflowNbr(
            dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            const CellID& cellID,
            const bool& copyMomentsOnly,
            const uint popID,
            const bool calculate_V_moments
         );
         void vlasovBoundaryCopyFromAllClosestNbrs(
            dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            const CellID& cellID,
            const uint popID,
            const bool calculate_V_moments
         );
         void vlasovBoundaryFluffyCopyFromAllCloseNbrs(
            dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
            const CellID& cellID,
            const uint popID,
            const bool calculate_V_moments,
            creal fluffiness
         );
         std::array<int, 3> getTheClosestNonsysboundaryCell(
            fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
            cint i,
            cint j,
            cint k
         );
         std::vector< std::array<int, 3> > getAllClosestNonsysboundaryCells(
            fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
            cint i,
            cint j,
            cint k
         );
         CellID & getTheClosestNonsysboundaryCell(
            const CellID& cellID
         );
         std::vector<CellID> & getAllClosestNonsysboundaryCells(
            const CellID& cellID
         );
         std::vector<CellID> & getAllCloseNonsysboundaryCells(
            const CellID& cellID
         );
         Real
         fieldBoundaryCopyFromSolvingNbrMagneticField(fsgrids::perbspan b,
                                                      fsgrids::consttechnicalspan technical,
                                                      const fsgrid::FsStencil& stencil, cuint component, cuint mask);
 
         CellID & getTheClosestL1OutflowCell(
            const CellID& cellID
         );
         std::vector<CellID> & getAllClosestL1OutflowCells(
            const CellID& cellID
         );

         uint precedence;
         bool dynamic;
         std::array<bool, 3> periodic;
         std::unordered_map<CellID, std::vector<CellID>> allClosestNonsysboundaryCells;
         std::unordered_map<CellID, std::vector<CellID>> allCloseNonsysboundaryCells;
         std::unordered_map<CellID, std::vector<CellID>> allClosestL1OutflowCells;
         std::unordered_map<CellID, std::vector<CellID>> allCloseL1OutflowCells;

         bool applyUponRestart;
   };
 
   class OuterBoundaryCondition: public SysBoundaryCondition {
      public:
         virtual void assignSysBoundary(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid, fsgrids::technicalspan technical, FieldSolverGrid &fsgrid);
      protected:
         std::array<bool, 6> facesToProcess;
   };
 
   // Moved outside the class since it's a helper function that doesn't require member access
   void averageCellData (
      const dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
      std::vector<CellID> cellList,
      SpatialCell *to,
      const uint popID,
      const creal fluffiness = 0
   );

   vmesh::LocalID findMaxwellianBlocksToInitialize(
      const uint popID,
      spatial_cell::SpatialCell& cell,
      creal& rho,
      creal& T,
      creal& VX0,
      creal& VY0,
      creal& VZ0);
 
   void determineFaceNoClassMembers(
      bool* isThisCellOnAFace,
      creal x, creal y, creal z,
      creal dx, creal dy, creal dz,
      const std::array<bool, 3> periodicity,
      const bool excludeSlicesAndPeriodicDimensions=false // (default)
   );
 
} // namespace SBC
 
#endif