sysboundarycondition.cpp

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.
 */

 
#include <cstdlib>
#include <iostream>
 
#include "../parameters.h"
#include "../vlasovsolver/vlasovmover.h"
#include "sysboundarycondition.h"
#include "../projects/projects_common.h"
#include "../object_wrapper.h"
 
using namespace std;
 
namespace SBC {
   // ************************************************************
   // ***** DEFINITIONS FOR BOUNDARYCONDITION BASE CLASS *****
   // ************************************************************

   void SysBoundaryCondition::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 (default)
   ) const {
      for(uint i=0; i<6; i++) {
         isThisCellOnAFace[i] = false;
      }
      if(x > Parameters::xmax - dx * 2) {
         isThisCellOnAFace[0] = true;
      }
      if(x < Parameters::xmin + dx * 2) {
         isThisCellOnAFace[1] = true;
      }
      if(y > Parameters::ymax - dy * 2) {
         isThisCellOnAFace[2] = true;
      }
      if(y < Parameters::ymin + dy * 2) {
         isThisCellOnAFace[3] = true;
      }
      if(z > Parameters::zmax - dz * 2) {
         isThisCellOnAFace[4] = true;
      }
      if(z < Parameters::zmin + dz * 2) {
         isThisCellOnAFace[5] = true;
      }
      if(excludeSlicesAndPeriodicDimensions == true) {
         if (Parameters::xcells_ini == 1 || this->periodic[0]) {
            isThisCellOnAFace[0] = false;
            isThisCellOnAFace[1] = false;
         }
         if (Parameters::ycells_ini == 1 || this->periodic[1]) {
            isThisCellOnAFace[2] = false;
            isThisCellOnAFace[3] = false;
         }
         if (Parameters::zcells_ini == 1 || this->periodic[2]) {
            isThisCellOnAFace[4] = false;
            isThisCellOnAFace[5] = false;
         }
      }
   }

   void SysBoundaryCondition::determineFace(
      std::array<bool, 6> &isThisCellOnAFace,
      const dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
      const CellID id,
      const bool excludeSlicesAndPeriodicDimensions //=false (default)
   ) {
      isThisCellOnAFace.fill(false);
 
      std::array<std::set<CellID>, 6> dirNeighbors;
      auto* p = mpiGrid.get_neighbors_of(id, Neighborhoods::VLASOV_SOLVER);
      if (!p) {
         std::cerr << "No neighbors found for " << id << std::endl;
         return;
      }
      int nbrs {0};
      for (auto pair : *p) {
         if (!pair.first) {
            continue;   // Error cells should obviously not be counted
         } else if (pair.second[0] > 0) {
            dirNeighbors[0].insert(pair.first);
         } else if (pair.second[0] < 0) {
            dirNeighbors[1].insert(pair.first);
         } else if (pair.second[1] > 0) {
            dirNeighbors[2].insert(pair.first);
         } else if (pair.second[1] < 0) {
            dirNeighbors[3].insert(pair.first);
         } else if (pair.second[2] > 0) {
            dirNeighbors[4].insert(pair.first);
         } else if (pair.second[2] < 0) {
            dirNeighbors[5].insert(pair.first);
         }
         ++nbrs;
         //std::cerr << pair.second[0] << " " << pair.second[1] << " " << pair.second[2] << std::endl;
      }
 
      for (int i = 0; i < 6; ++i) {
         if(dirNeighbors[i].size() < 2) {
            isThisCellOnAFace[i] = true;
         }
      }
 
      if(excludeSlicesAndPeriodicDimensions == true) {
         if(Parameters::xcells_ini == 1 || this->periodic[0]) {
            isThisCellOnAFace[0] = false;
            isThisCellOnAFace[1] = false;
         }
         if(Parameters::ycells_ini == 1 || this->periodic[1]) {
            isThisCellOnAFace[2] = false;
            isThisCellOnAFace[3] = false;
         }
         if(Parameters::zcells_ini == 1 || this->periodic[2]) {
            isThisCellOnAFace[4] = false;
            isThisCellOnAFace[5] = false;
         }
      }
      return;
   }

   SysBoundaryCondition::SysBoundaryCondition() { }

   SysBoundaryCondition::~SysBoundaryCondition() { }

   void SysBoundaryCondition::addParameters() {
      cerr << "ERROR: SysBoundaryCondition::addParameters called instead of derived class function!" << endl;
   }

   void
   SysBoundaryCondition::setCellDerivativesToZero(fsgrids::dperbspan dperb,
                                                  fsgrids::dmomentsspan dmoments,
                                                  const fsgrid::FsStencil& stencil, cuint component) {
      auto& dPerBGrid0 = dperb[stencil.ooo()];
      auto& dMomentsGrid0 = dmoments[stencil.ooo()];
      switch(component) {
         case 0: // x, xx
            dMomentsGrid0[fsgrids::dmoments::drhomdx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::drhoqdx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp11dx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp22dx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp33dx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVxdx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVydx] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVzdx] = 0.0;
 
            dPerBGrid0[fsgrids::dperb::dPERBydx] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBzdx] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBydxx] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBzdxx] = 0.0;
            break;
         case 1: // y, yy
            dMomentsGrid0[fsgrids::dmoments::drhomdy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::drhoqdy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp11dy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp22dy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp33dy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVxdy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVydy] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVzdy] = 0.0;
 
            dPerBGrid0[fsgrids::dperb::dPERBxdy] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBzdy] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBxdyy] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBzdyy] = 0.0;
            break;
         case 2: // z, zz
            dMomentsGrid0[fsgrids::dmoments::drhomdz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::drhoqdz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp11dz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp22dz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dp33dz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVxdz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVydz] = 0.0;
            dMomentsGrid0[fsgrids::dmoments::dVzdz] = 0.0;
 
            dPerBGrid0[fsgrids::dperb::dPERBxdz] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBydz] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBxdzz] = 0.0;
            dPerBGrid0[fsgrids::dperb::dPERBydzz] = 0.0;
            break;
         case 3: // xy
            dPerBGrid0[fsgrids::dperb::dPERBzdxy] = 0.0;
            break;
         case 4: // xz
            dPerBGrid0[fsgrids::dperb::dPERBydxz] = 0.0;
            break;
         case 5: // yz
            dPerBGrid0[fsgrids::dperb::dPERBxdyz] = 0.0;
            break;
         default:
            cerr << __FILE__ << ":" << __LINE__ << ":" << " Invalid component" << endl;
            abort_mpi("Invalid component", 1);
      }
   }

   void SysBoundaryCondition::setCellBVOLDerivativesToZero(fsgrids::volspan vols,
                                                           const fsgrid::FsStencil& stencil, cuint component) {
      auto& vol = vols[stencil.ooo()];
      switch(component) {
         case 0:
            vol[fsgrids::volfields::dPERBXVOLdx] = 0.0;
            vol[fsgrids::volfields::dPERBYVOLdx] = 0.0;
            vol[fsgrids::volfields::dPERBZVOLdx] = 0.0;
            break;
         case 1:
            vol[fsgrids::volfields::dPERBXVOLdy] = 0.0;
            vol[fsgrids::volfields::dPERBYVOLdy] = 0.0;
            vol[fsgrids::volfields::dPERBZVOLdy] = 0.0;
            break;
         case 2:
            vol[fsgrids::volfields::dPERBXVOLdz] = 0.0;
            vol[fsgrids::volfields::dPERBYVOLdz] = 0.0;
            vol[fsgrids::volfields::dPERBZVOLdz] = 0.0;
            break;
         default:
         abort_mpi("Invalid component", 1);
      }
   }

   void SysBoundaryCondition::vlasovBoundaryCopyFromTheClosestNbr(
         dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
         const CellID& cellID,
         const bool& copyMomentsOnly,
         const uint popID,
         const bool copy_V_moments
   ) {
      const CellID closestCell = getTheClosestNonsysboundaryCell(cellID);
 
      if(closestCell == INVALID_CELLID) {
         abort_mpi("No closest cell found!", 1);
      }
      copyCellData(mpiGrid[closestCell],mpiGrid[cellID], copyMomentsOnly, popID, copy_V_moments);
   }

   void SysBoundaryCondition::vlasovBoundaryCopyFromTheClosestL1OutflowNbr(
         dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
         const CellID& cellID,
         const bool& copyMomentsOnly,
         const uint popID,
         const bool copy_V_moments
   ) {
      const CellID closestCell = getTheClosestL1OutflowCell(cellID);
 
      if(closestCell == INVALID_CELLID) {
         abort_mpi("No closest L1 Outflow cell found!", 1);
      }
      phiprof::Timer boundaryTimer {"vlasovBoundaryCopyFromTheClosestL1OutflowNbr"};
      copyCellData(mpiGrid[closestCell],mpiGrid[cellID], copyMomentsOnly, popID, copy_V_moments);
      boundaryTimer.stop();
   }

   void SysBoundaryCondition::vlasovBoundaryCopyFromAllClosestNbrs(
      dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
      const CellID& cellID,const uint popID, const bool copy_V_moments
   ) {
      const vector<CellID>& closestCells = getAllClosestNonsysboundaryCells(cellID);
 
      if(closestCells[0] == INVALID_CELLID) {
         abort_mpi("No closest cell found!", 1);
      }
      averageCellData(mpiGrid, closestCells, mpiGrid[cellID], popID);
   }

   void SysBoundaryCondition::vlasovBoundaryFluffyCopyFromAllCloseNbrs(
      dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
      const CellID& cellID,const uint popID,const bool copy_V_moments,const creal fluffiness
   ) {
      const vector<CellID>& closeCells = getAllCloseNonsysboundaryCells(cellID);
 
      if(closeCells[0] == INVALID_CELLID) {
         abort_mpi("No close cell found!", 1);
      }
      averageCellData(mpiGrid, closeCells, mpiGrid[cellID], popID, fluffiness);
   }

   void SysBoundaryCondition::copyCellData(
            const SpatialCell* const from,
            SpatialCell* to,
            const bool copyMomentsOnly,
            const uint popID,
            const bool copy_V_moments
   ) {
      if (popID == 0) {
         if (copy_V_moments) {
            to->parameters[CellParams::RHOM_V] = from->parameters[CellParams::RHOM_V];
            to->parameters[CellParams::VX_V] = from->parameters[CellParams::VX_V];
            to->parameters[CellParams::VY_V] = from->parameters[CellParams::VY_V];
            to->parameters[CellParams::VZ_V] = from->parameters[CellParams::VZ_V];
            to->parameters[CellParams::RHOQ_V] = from->parameters[CellParams::RHOQ_V];
            to->parameters[CellParams::P_11_V] = from->parameters[CellParams::P_11_V];
            to->parameters[CellParams::P_22_V] = from->parameters[CellParams::P_22_V];
            to->parameters[CellParams::P_33_V] = from->parameters[CellParams::P_33_V];
         } else {
            to->parameters[CellParams::RHOM_R] = from->parameters[CellParams::RHOM_R];
            to->parameters[CellParams::VX_R] = from->parameters[CellParams::VX_R];
            to->parameters[CellParams::VY_R] = from->parameters[CellParams::VY_R];
            to->parameters[CellParams::VZ_R] = from->parameters[CellParams::VZ_R];
            to->parameters[CellParams::RHOQ_R] = from->parameters[CellParams::RHOQ_R];
            to->parameters[CellParams::P_11_R] = from->parameters[CellParams::P_11_R];
            to->parameters[CellParams::P_22_R] = from->parameters[CellParams::P_22_R];
            to->parameters[CellParams::P_33_R] = from->parameters[CellParams::P_33_R];
         }
      }
 
      if(copyMomentsOnly) {
         to->get_population(popID).RHO = from->get_population(popID).RHO;
         if (copy_V_moments) {
            to->get_population(popID).RHO_V = from->get_population(popID).RHO_V;
         } else {
            to->get_population(popID).RHO_R = from->get_population(popID).RHO_R;
         }
         for (uint i=0; i<3; i++) {
            to->get_population(popID).V[i] = from->get_population(popID).V[i];
            if (copy_V_moments) {
               to->get_population(popID).V_V[i] = from->get_population(popID).V_V[i];
            } else {
               to->get_population(popID).V_R[i] = from->get_population(popID).V_R[i];
            }
         }
         for (uint i=0; i<6; i++) {
            to->get_population(popID).P[i] = from->get_population(popID).P[i];
            if (copy_V_moments) {
               to->get_population(popID).P_V[i] = from->get_population(popID).P_V[i];
            } else {
               to->get_population(popID).P_R[i] = from->get_population(popID).P_R[i];
            }
         }
      } else {
         to->set_population(from->get_population(popID), popID);
      }
   }

   void averageCellData(
         const dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
         const vector<CellID> cellList,
         SpatialCell *to,
         const uint popID,
         const creal fluffiness /* default =0.0*/
   ) {
      const size_t numberOfCells = cellList.size();
      creal factor = fluffiness / convert<Real>(numberOfCells);
 
      // Rescale own vspace
      if (fluffiness != 0.0) {
         to->scale_population(1.0 - fluffiness, popID);
      }
      // Add other cell vspaces
      for (size_t i=0; i<numberOfCells; i++) {
         const SpatialCell* from = mpiGrid[cellList[i]];
         to->increment_population(from->get_population(popID), factor, popID);
      }
   }

   bool SysBoundaryCondition::updateSysBoundaryConditionsAfterLoadBalance(
      const dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
      const vector<CellID> & local_cells_on_boundary
   ) {
      // Loop over cellids
      for(vector<CellID>::const_iterator it = local_cells_on_boundary.begin(); it != local_cells_on_boundary.end(); ++it) {
         const CellID cellId = *it;
         vector<CellID> & closestCells = allClosestNonsysboundaryCells[cellId];
         closestCells.clear();
         vector<CellID> & closeCells = allCloseNonsysboundaryCells[cellId];
         closeCells.clear();
         vector<CellID> & closestL1OutflowCells = allClosestL1OutflowCells[cellId];
         closestL1OutflowCells.clear();
         vector<CellID> & closeL1OutflowCells = allCloseL1OutflowCells[cellId];
         closeL1OutflowCells.clear();
         uint dist = numeric_limits<uint>::max();
         uint d2 = numeric_limits<uint>::max();
         uint distL1 = numeric_limits<uint>::max();
         uint d2L1 = numeric_limits<uint>::max();
 
         // Only closer neighborhood for layer 1
         if(mpiGrid[cellId]->sysBoundaryLayer == 1) {
            for (auto nbrPair : *mpiGrid.get_neighbors_of(cellId, Neighborhoods::SYSBOUNDARIES)) {
               if(nbrPair.first != INVALID_CELLID) {
                  CellID neighbor = nbrPair.first;
                  if(mpiGrid[neighbor]->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY) {
                     // Find distance and update closestCells
                     d2 = nbrPair.second[0]*nbrPair.second[0] + nbrPair.second[1]*nbrPair.second[1] + nbrPair.second[2]*nbrPair.second[2];
                     for (auto i : *mpiGrid.get_neighbors_to(cellId, Neighborhoods::SYSBOUNDARIES)) {
                        if (i.first == neighbor) {
                           closeCells.push_back(neighbor);
                           if(d2 == dist) {
                              closestCells.push_back(neighbor);
                           } else if (d2 < dist) {
                              closestCells.clear();
                              closestCells.push_back(neighbor);
                              dist = d2;
                           }
                        }
                     }
                  }
               }
            }
         }
         // search further for L2
         if (mpiGrid[cellId]->sysBoundaryLayer == 2) {
            for (auto nbrPair : *mpiGrid.get_neighbors_of(cellId, Neighborhoods::SYSBOUNDARIES_EXTENDED)) {
               if(nbrPair.first != INVALID_CELLID) {
                  CellID neighbor = nbrPair.first;
                  if(mpiGrid[neighbor]->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY) {
                     // Find distance and update closestCells
                     d2 = nbrPair.second[0]*nbrPair.second[0] + nbrPair.second[1]*nbrPair.second[1] + nbrPair.second[2]*nbrPair.second[2];
                     closeCells.push_back(neighbor);
                     if(d2 == dist) {
                        closestCells.push_back(neighbor);
                     } else if (d2 < dist) {
                        closestCells.clear();
                        closestCells.push_back(neighbor);
                        dist = d2;
                     }
                  }
                  if(mpiGrid[cellId]->sysBoundaryFlag == sysboundarytype::OUTFLOW && mpiGrid[neighbor]->sysBoundaryFlag == sysboundarytype::OUTFLOW && mpiGrid[neighbor]->sysBoundaryLayer == 1) {
                     // Find distance and update closestCells
                     d2L1 = nbrPair.second[0]*nbrPair.second[0] + nbrPair.second[1]*nbrPair.second[1] + nbrPair.second[2]*nbrPair.second[2];
                     closeL1OutflowCells.push_back(neighbor);
                     if(d2L1 == distL1) {
                        closestL1OutflowCells.push_back(neighbor);
                     } else if (d2L1 < distL1) {
                        closestL1OutflowCells.clear();
                        closestL1OutflowCells.push_back(neighbor);
                        distL1 = d2L1;
                     }
                  }
               }
            }
         }
 
         if(closestCells.size() == 0) {
            closestCells.push_back(INVALID_CELLID);
         }
         if(closeCells.size() == 0) {
            closeCells.push_back(INVALID_CELLID);
         }
         if(closestL1OutflowCells.size() == 0) {
            closestL1OutflowCells.push_back(INVALID_CELLID);
         }
         if(closeL1OutflowCells.size() == 0) {
            closeL1OutflowCells.push_back(INVALID_CELLID);
         }
      }
      return true;
   }

   array<int, 3> SysBoundaryCondition::getTheClosestNonsysboundaryCell(
      fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
      cint i,
      cint j,
      cint k
   ) {
      const vector< array<int, 3> > closestCells = getAllClosestNonsysboundaryCells(technical, fsgrid, i, j, k);
      return closestCells.at(0);
   }

   vector< array<int, 3> > SysBoundaryCondition::getAllClosestNonsysboundaryCells(
      fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
      cint i,
      cint j,
      cint k
   ) {
      const auto stencil = fsgrid.makeStencil(i, j, k);
      int distance = numeric_limits<int>::max();
      vector<array<int, 3>> closestCells;
 
      for (int kk=-2; kk<3; kk++) {
         for (int jj=-2; jj<3; jj++) {
            for (int ii=-2; ii<3 ; ii++) {
               if (stencil.cellExists(ii, jj, kk)) {
                  if (technical[stencil.indexFromOffset(ii, jj, kk)].sysBoundaryFlag ==
                      sysboundarytype::NOT_SYSBOUNDARY) {
                     distance = min(distance, ii * ii + jj * jj + kk * kk);
                  }
               }
            }
         }
      }
 
      for (int kk=-2; kk<3; kk++) {
         for (int jj=-2; jj<3; jj++) {
            for (int ii=-2; ii<3 ; ii++) {
               if (stencil.cellExists(ii, jj, kk)) {
                  if (technical[stencil.indexFromOffset(ii, jj, kk)].sysBoundaryFlag ==
                      sysboundarytype::NOT_SYSBOUNDARY) {
                     const int d = ii * ii + jj * jj + kk * kk;
                     if (d == distance) {
                        closestCells.push_back({i + ii, j + jj, k + kk});
                     }
                  }
               }
            }
         }
      }
 
      if (closestCells.empty()) {
         closestCells.push_back({numeric_limits<int>::min()});
      }
 
      return closestCells;
   }

   CellID & SysBoundaryCondition::getTheClosestNonsysboundaryCell(
      const CellID& cellID
   ) {
      vector<CellID> & closestCells = allClosestNonsysboundaryCells.at(cellID);
      return closestCells.at(0);
   }

   CellID & SysBoundaryCondition::getTheClosestL1OutflowCell(
      const CellID& cellID
   ) {
      vector<CellID> & closestCells = allClosestL1OutflowCells.at(cellID);
      return closestCells.at(0);
   }

   vector<CellID> & SysBoundaryCondition::getAllClosestNonsysboundaryCells(
      const CellID& cellID
   ) {
      vector<CellID> & closestCells = allClosestNonsysboundaryCells.at(cellID);
      return closestCells;
   }

   vector<CellID> & SysBoundaryCondition::getAllCloseNonsysboundaryCells(
      const CellID& cellID
   ) {
      vector<CellID> & closeCells = allCloseNonsysboundaryCells.at(cellID);
      return closeCells;
   }
 
   Real SysBoundaryCondition::fieldBoundaryCopyFromSolvingNbrMagneticField(
       fsgrids::perbspan b, fsgrids::consttechnicalspan technical,
       const fsgrid::FsStencil& stencil, cuint component, cuint mask) {
      int distance = numeric_limits<int>::max();
      auto closestCellIndex = 0;
 
      for (auto kk = -2; kk < 3; kk++) {
         for (auto jj = -2; jj < 3; jj++) {
            for (auto ii = -2; ii < 3; ii++) {
               if (stencil.cellExists(ii, jj, kk)) {
                  const auto index = stencil.indexFromOffset(ii, jj, kk);
                  const auto& tech = technical[index];
                  const bool copyable = (tech.SOLVE & mask) == mask &&
                                        tech.sysBoundaryFlag != sysboundarytype::DO_NOT_COMPUTE &&
                                        tech.sysBoundaryFlag != sysboundarytype::OUTER_BOUNDARY_PADDING;
                  const int d = ii * ii + jj * jj + kk * kk;
                  if (copyable && d < distance) {
                     distance = d;
                     closestCellIndex = index;
                  }
               }
            }
         }
      }
 
      if (distance == numeric_limits<int>::max()) {
         abort_mpi("No closest cell found!", 1);
      }
 
      return b[closestCellIndex][fsgrids::bfield::PERBX + component];
   }

   void SysBoundaryCondition::getFaces(bool* faces) {
      cerr << "ERROR: SysBoundaryCondition::getFaces called instead of derived class function!" << endl;
      for(int i=0; i<6; i++) {
        faces[i]=false;
      }
   }

   uint SysBoundaryCondition::getPrecedence() const {return precedence;}

   bool SysBoundaryCondition::isDynamic() const { return dynamic; }
   
   void SysBoundaryCondition::setPeriodicity(
      std::array<bool, 3> isFacePeriodic
   ) {
      this->periodic = isFacePeriodic;
   }
   
   bool SysBoundaryCondition::doApplyUponRestart() const {return this->applyUponRestart;}
   
   void SysBoundaryCondition::mapCellPotentialAndGetEXBDrift(
      std::array<Real, CellParams::N_SPATIAL_CELL_PARAMS>& cellParams
   ) {
      std::cerr << "Error: SysBoundaryCondition::mapCellPotentialAndGetEXBDrift called!\n";
   }
 
   vmesh::LocalID findMaxwellianBlocksToInitialize(
      const uint popID,
      spatial_cell::SpatialCell& cell,
      creal& rho,
      creal& T,
      creal& VX0,
      creal& VY0,
      creal& VZ0) {
 
      bool search = true;
      uint counter = 0;
      const vmesh::LocalID* vblocks_ini = cell.get_velocity_grid_length(popID);
      vmesh::VelocityMesh *vmesh = cell.get_velocity_mesh(popID);
      const Real mass = getObjectWrapper().particleSpecies[popID].mass;
 
      vmesh::GlobalID *GIDbuffer;
      #ifdef USE_GPU
      // Host-pinned memory buffer, max possible size
      const uint blocksCount = vblocks_ini[0]*vblocks_ini[1]*vblocks_ini[2];
      CHK_ERR( gpuMallocHost((void**)&GIDbuffer,blocksCount*sizeof(vmesh::GlobalID)) );
      #endif
      // Non-GPU: insert directly into vmesh
 
      Real V_crds[3];
      Real dV[3];
      dV[0] = cell.get_velocity_grid_block_size(popID)[0];
      dV[1] = cell.get_velocity_grid_block_size(popID)[1];
      dV[2] = cell.get_velocity_grid_block_size(popID)[2];
      creal minValue = cell.getVelocityBlockMinValue(popID);
      // Single cell, not block
      const Real dvx=cell.get_velocity_grid_cell_size(popID)[0];
      const Real dvy=cell.get_velocity_grid_cell_size(popID)[1];
      const Real dvz=cell.get_velocity_grid_cell_size(popID)[2];
 
      while (search) {
         if (0.1 * minValue > projects::MaxwellianPhaseSpaceDensity(counter*dV[0]+0.5*dvx, 0.5*dvy, 0.5*dvz, T, rho, mass) || counter > vblocks_ini[0]) {
            search = false;
         }
         counter++;
      }
      counter+=2;
 
      Real vRadiusSquared = (Real)counter * (Real)counter * dV[0] * dV[0];
 
      #ifndef USE_GPU
      // sphere volume is 4/3 pi r^3, approximate that 5*counterX*counterY*counterZ is enough.
      vmesh::LocalID currentMaxSize = 5*counter*counter*counter;
      vmesh->setNewSize(currentMaxSize);
      GIDbuffer = vmesh->getGrid()->data();
      #endif
 
      vmesh::LocalID LID = 0;
      for (uint kv=0; kv<vblocks_ini[2]; ++kv) {
         for (uint jv=0; jv<vblocks_ini[1]; ++jv) {
            for (uint iv=0; iv<vblocks_ini[0]; ++iv) {
               const vmesh::GlobalID GID = vmesh->getGlobalID(iv,jv,kv);
 
               cell.get_velocity_block_coordinates(popID,GID,V_crds);
               V_crds[0] += 0.5*dV[0] - VX0;
               V_crds[1] += 0.5*dV[1] - VY0;
               V_crds[2] += 0.5*dV[2] - VZ0;
               Real R2 = ((V_crds[0])*(V_crds[0])
                          + (V_crds[1])*(V_crds[1])
                          + (V_crds[2])*(V_crds[2]));
 
               #ifndef USE_GPU
               if (LID >= currentMaxSize) {
                  currentMaxSize = LID + counter*counter*counter;
                  vmesh->setNewSize(currentMaxSize);
                  GIDbuffer = vmesh->getGrid()->data();
               }
               #endif
               if (R2 < vRadiusSquared) {
                  GIDbuffer[LID] = GID;
                  LID++;
               }
            }
         }
      }
 
      // Set final size of vmesh
      cell.get_population(popID).N_blocks = LID;
 
      #ifdef USE_GPU
      // Copy data into place
      cell.dev_resize_vmesh(popID,LID);
      vmesh::GlobalID *GIDtarget = vmesh->getGrid()->data();
      gpuStream_t stream = gpu_getStream();
      CHK_ERR( gpuMemcpyAsync(GIDtarget, GIDbuffer, LID*sizeof(vmesh::GlobalID), gpuMemcpyHostToDevice, stream));
      CHK_ERR( gpuStreamSynchronize(stream) );
      CHK_ERR( gpuFreeHost(GIDbuffer));
      #else
      vmesh->setNewSize(LID);
      #endif
 
      return LID;
   }
 
   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)
   ) {
      for(uint i=0; i<6; i++) {
         isThisCellOnAFace[i] = false;
      }
      if(x > Parameters::xmax - dx * 2) {
         isThisCellOnAFace[0] = true;
      }
      if(x < Parameters::xmin + dx * 2) {
         isThisCellOnAFace[1] = true;
      }
      if(y > Parameters::ymax - dy * 2) {
         isThisCellOnAFace[2] = true;
      }
      if(y < Parameters::ymin + dy * 2) {
         isThisCellOnAFace[3] = true;
      }
      if(z > Parameters::zmax - dz * 2) {
         isThisCellOnAFace[4] = true;
      }
      if(z < Parameters::zmin + dz * 2) {
         isThisCellOnAFace[5] = true;
      }
      if(excludeSlicesAndPeriodicDimensions == true) {
         if (Parameters::xcells_ini == 1 || periodicity[0]) {
            isThisCellOnAFace[0] = false;
            isThisCellOnAFace[1] = false;
         }
         if (Parameters::ycells_ini == 1 || periodicity[1]) {
            isThisCellOnAFace[2] = false;
            isThisCellOnAFace[3] = false;
         }
         if (Parameters::zcells_ini == 1 || periodicity[2]) {
            isThisCellOnAFace[4] = false;
            isThisCellOnAFace[5] = false;
         }
      }
   }
 
   // The function above and the commented parts below make an implementation with lambda
   // but it breaks some communicator and YPK did not manage to get a working version out of this.
   // Left in the copies of variable so the [=] compiles without complaining.
 
   void OuterBoundaryCondition::assignSysBoundary(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                  fsgrids::technicalspan technical, FieldSolverGrid &fsgrid) {
      bool doAssign;
      std::array<bool, 6> isThisCellOnAFace;
 
      // Assign boundary flags to local DCCRG cells
      const vector<CellID>& cells = getLocalCells();
      for (auto& id : cells) {
         if (mpiGrid[id]->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
            continue;
         }
         creal* const cellParams = &(mpiGrid[id]->parameters[0]);
         creal dx = cellParams[CellParams::DX];
         creal dy = cellParams[CellParams::DY];
         creal dz = cellParams[CellParams::DZ];
         creal x = cellParams[CellParams::XCRD] + 0.5 * dx;
         creal y = cellParams[CellParams::YCRD] + 0.5 * dy;
         creal z = cellParams[CellParams::ZCRD] + 0.5 * dz;
 
         for(uint i=0; i<6; i++) {
            isThisCellOnAFace[i] = false;
         }
         determineFace(isThisCellOnAFace.data(), x, y, z, dx, dy, dz);
         // Comparison of the array defining which faces to use and the array telling on which faces this cell is
         doAssign = false;
         for (int j = 0; j < 6; j++) {
            doAssign = doAssign || (facesToProcess[j] && isThisCellOnAFace[j]);
         }
         if (doAssign) {
            mpiGrid[id]->sysBoundaryFlag = this->getIndex();
         }
      }
      const auto index_local = this->getIndex();
      const std::array<bool, 6> facesToProcess_local = facesToProcess;
      const std::array<bool, 3> periodic_local = this->periodic;
      // Assign boundary flags to local fsgrid cells
      fsgrid.serial_for([](int timerId) -> phiprof::Timer { return phiprof::Timer{timerId}; },
                        phiprof::initializeTimer("Assign sysboundary flags to fsgrid cells"), technical,
                        [=](const fsgrid::Coordinates &coordinates, const fsgrid::FsStencil& stencil, cuint sysBoundaryFlag, cuint sysBoundaryLayer) {
         creal dx = P::dx_ini * pow(2, -technical[stencil.ooo()].refLevel);
         creal dy = P::dy_ini * pow(2, -technical[stencil.ooo()].refLevel);
         creal dz = P::dz_ini * pow(2, -technical[stencil.ooo()].refLevel);
 
         std::array<bool, 6> isThisCellOnAFace = {{false}};
         bool doAssign = false;
         const std::array<Real, 3> coords = coordinates.getPhysicalCoords(stencil.i, stencil.j, stencil.k);
         const std::array<Real, 3> gridSpacing = coordinates.physicalGridSpacing;
 
         determineFaceNoClassMembers(isThisCellOnAFace.data(), coords[0] + 0.5 * gridSpacing[0], coords[1] + 0.5 * gridSpacing[1], coords[2] + 0.5 * gridSpacing[2], dx, dy, dz, periodic_local);
         for (int iface = 0; iface < 6; iface++) {
            doAssign = doAssign || (facesToProcess_local[iface] && isThisCellOnAFace[iface]);
         }
         if (doAssign) {
            technical[stencil.ooo()].sysBoundaryFlag = index_local;
         }
      });
   }
 
} // namespace SBC