SysBoundary communicate VDFs only in smaller neighborhood

datareduction/datareducer.cpp

@@ -314,6 +314,10 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
             const std::string& pop = species.name;
             outputReducer->addOperator(new DRO::DataReductionOperatorPopulations<Real>(pop + "/vg_rho", i, offsetof(spatial_cell::Population, RHO), 1));
             outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$n_\\mathrm{"+pop+"}$","1.0");
+            outputReducer->addOperator(new DRO::DataReductionOperatorPopulations<Real>(pop + "/vg_rho_r", i, offsetof(spatial_cell::Population, RHO_R), 1));
+            outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$n_\\mathrm{"+pop+"}$","1.0");
+            outputReducer->addOperator(new DRO::DataReductionOperatorPopulations<Real>(pop + "/vg_rho_v", i, offsetof(spatial_cell::Population, RHO_V), 1));
+            outputReducer->addMetadata(outputReducer->size()-1,"1/m^3","$\\mathrm{m}^{-3}$","$n_\\mathrm{"+pop+"}$","1.0");
          }
          if(!P::systemWriteAllDROs) {
             continue;

grid.cpp

@@ -364,7 +364,12 @@ void initializeGrids(
       calculateInitialVelocityMoments(mpiGrid);
    } else {
       phiprof::Timer timer {"Init moments"};
+      #pragma omp parallel for schedule(guided,1)
       for (size_t i=0; i<cells.size(); ++i) {
+         // easier to skip here than adding one more bool flag to calculateCellMoments - handles L2 outflow cells without VDF
+         if(mpiGrid[cells[i]]->sysBoundaryFlag == sysboundarytype::OUTFLOW && mpiGrid[cells[i]]->sysBoundaryLayer != 1) {
+            continue;
+         }
          calculateCellMoments(mpiGrid[cells[i]], true, true);
       }
    }

sysboundary/copysphere.cpp

@@ -832,6 +832,20 @@ namespace SBC {
 
       calculateCellMoments(&templateCell,true,false,true);
 
+      for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
+         Population &pop = templateCell.get_population(popID);
+         pop.RHO_V = pop.RHO;
+         pop.RHO_R = pop.RHO;
+         for(int i=0; i<3; i++) {
+            pop.V_R[i] = pop.V[i];
+            pop.V_V[i] = pop.V[i];
+         }
+         for(int i=0; i<6; i++) {
+            pop.P_R[i] = pop.P[i];
+            pop.P_V[i] = pop.P[i];
+         }
+      }
+
       // WARNING Time-independence assumed here. Normal moments computed in setProjectCell
       templateCell.parameters[CellParams::RHOM_R] = templateCell.parameters[CellParams::RHOM];
       templateCell.parameters[CellParams::VX_R] = templateCell.parameters[CellParams::VX];

sysboundary/copysphere.h

@@ -132,6 +132,9 @@ namespace SBC {
          const uint popID,
          const bool calculate_V_moments
       ) override;
+      virtual void setupL2OutflowAtRestart(
+         dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid
+      ) override {std::cerr << "ERROR: Copysphere::setupL2OutflowAtRestart called!" << std::endl;}
 
       void getFaces(bool *faces) override;
       virtual std::string getName() const override;

sysboundary/donotcompute.h

@@ -121,6 +121,10 @@ namespace SBC {
           const uint popID,
           const bool calculate_V_moments
       ) override { std::cerr << "ERROR: DoNotCompute::vlasovBoundaryCondition called!" << std::endl;}
+      virtual void setupL2OutflowAtRestart(
+         dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid
+      ) override {std::cerr << "ERROR: DoNotCompute::setupL2OutflowAtRestart called!" << std::endl;}
+
    };
 }
 

sysboundary/inflow.h

@@ -94,6 +94,9 @@ class Inflow : public OuterBoundaryCondition {
                                           cint i, cint j, cint k, cuint component) override;
    virtual void vlasovBoundaryCondition(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
                                         const CellID& cellID, const uint popID, const bool doCalcMomentsV) override;
+   virtual void setupL2OutflowAtRestart(
+      dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid
+   ) override {std::cerr << "ERROR: Inflow::setupL2OutflowAtRestart called!" << std::endl;}
 
    virtual void getFaces(bool* faces) override;
    virtual std::string getName() const override = 0;

sysboundary/ionosphere.cpp

@@ -3381,8 +3381,8 @@ namespace SBC {
          // let's get rid of blocks not fulfilling the criteria here to save memory.
          mpiGrid[cellID]->adjustSingleCellVelocityBlocks(popID,true);
 
-         // TODO: The moments can also be analytically calculated from ionosphere parameters.
-         // Maybe that's faster?
+         // In principle this could call _R or _V instead according to calculate_V_moments (unused at the moment)
+         // But the relevant moments will get recomputed in other spots when needed.
          calculateCellMoments(mpiGrid[cellID], true, false, true);
       } // End of if for coupling interval, we skip this altogether
 
@@ -3474,6 +3474,20 @@ namespace SBC {
 
       calculateCellMoments(&templateCell,true,false,true);
 
+      for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
+         Population &pop = templateCell.get_population(popID);
+         pop.RHO_V = pop.RHO;
+         pop.RHO_R = pop.RHO;
+         for(int i=0; i<3; i++) {
+            pop.V_R[i] = pop.V[i];
+            pop.V_V[i] = pop.V[i];
+         }
+         for(int i=0; i<6; i++) {
+            pop.P_R[i] = pop.P[i];
+            pop.P_V[i] = pop.P[i];
+         }
+      }
+
       // WARNING Time-independence assumed here. Normal moments computed in setProjectCell
       templateCell.parameters[CellParams::RHOM_R] = templateCell.parameters[CellParams::RHOM];
       templateCell.parameters[CellParams::VX_R] = templateCell.parameters[CellParams::VX];

sysboundary/ionosphere.h

@@ -54,7 +54,6 @@ namespace SBC {
       Real rho;
       Real V0[3];
       Real T;
-      Real fluffiness;
    };
 
    enum IonosphereBoundaryVDFmode { // How are inner boundary VDFs constructed from the ionosphere
@@ -395,6 +394,10 @@ namespace SBC {
          const uint popID,
          const bool calculate_V_moments
       ) override;
+      virtual void setupL2OutflowAtRestart(
+         dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid
+      ) override {std::cerr << "ERROR: Ionosphere::setupL2OutflowAtRestart called!" << std::endl;}
+
       virtual void updateState(
          dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid,
          FsGrid<std::array<Real, fsgrids::bfield::N_BFIELD>, FS_STENCIL_WIDTH>& perBGrid,

sysboundary/outflow.cpp

@@ -271,7 +271,7 @@ namespace SBC {
          }
 
          if (doApply) {
-            project.setCell(cell);
+            project.setCell(cell); // We set everything including VDF even in L2 cells to avoid a pile of spaghetti. Won't get communicated.
             cell->parameters[CellParams::RHOM_DT2] = cell->parameters[CellParams::RHOM];
             cell->parameters[CellParams::RHOQ_DT2] = cell->parameters[CellParams::RHOQ];
             cell->parameters[CellParams::VX_DT2] = cell->parameters[CellParams::VX];
@@ -411,7 +411,11 @@ namespace SBC {
                case vlasovscheme::NONE:
                   break;
                case vlasovscheme::COPY:
-                  vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,false,popID,calculate_V_moments);
+                  if(mpiGrid[cellID]->sysBoundaryLayer == 1) {
+                     vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,false,popID,calculate_V_moments); // copies VDF too
+                  } else {
+                     vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,true,popID,calculate_V_moments); // no VDF copy
+                  }
                   break;
                default:
                   abort_mpi("ERROR: invalid Outflow Vlasov scheme", 1);
@@ -420,6 +424,87 @@ namespace SBC {
       }
    }
 
+   void Outflow::setupL2OutflowAtRestart(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid) {
+      SpatialCell::set_mpi_transfer_type(Transfer::ALL_DATA);
+      mpiGrid.update_copies_of_remote_neighbors(Neighborhoods::SYSBOUNDARIES);
+
+      const vector<CellID>& cells = getLocalCells();
+      #pragma omp parallel
+      {
+         #pragma omp for schedule(guided,1)
+         for(uint i=0; i<cells.size(); i++) {
+            const CellID cellID = cells[i];
+            if(mpiGrid[cellID]->sysBoundaryFlag != this->getIndex()) {
+               continue;
+            }
+
+            std::array<bool, 6> isThisCellOnAFace;
+            determineFace(isThisCellOnAFace, mpiGrid, cellID, true);
+
+            for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
+               const OutflowSpeciesParameters& sP = this->speciesParams[popID];
+               for(uint i=0; i<6; i++) {
+                  if(isThisCellOnAFace[i] && facesToProcess[i] && !sP.facesToSkipVlasov[i]) {
+                     switch(sP.faceVlasovScheme[i]) {
+                        case vlasovscheme::NONE:
+                           break;
+                        case vlasovscheme::COPY:
+                           if(mpiGrid[cellID]->sysBoundaryLayer == 1) {
+                              vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,false,popID,true); // first false means copy VDF too, second true means V moments
+                              vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,false,popID,false); // first false means copy VDF too, second false means R moments
+                           }
+                           break;
+                        default:
+                           abort_mpi("ERROR: invalid Outflow Vlasov scheme", 1);
+                     } // switch
+                  } // if on face
+               } // faces
+            } // populations
+         } // cells
+
+         #pragma omp barrier
+         #pragma omp single
+         {
+            SpatialCell::set_mpi_transfer_type(Transfer::ALL_DATA);
+            mpiGrid.update_copies_of_remote_neighbors(Neighborhoods::SYSBOUNDARIES);
+         }
+         #pragma omp barrier
+
+         // then 2nd pass and copy from the closest L1 outflow neighbor
+         #pragma omp for schedule(guided,1)
+         for(uint i=0; i<cells.size(); i++) {
+            const CellID cellID = cells[i];
+            if(mpiGrid[cellID]->sysBoundaryFlag != this->getIndex()) {
+               continue;
+            }
+
+            std::array<bool, 6> isThisCellOnAFace;
+            determineFace(isThisCellOnAFace, mpiGrid, cellID, true);
+
+            for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
+               const OutflowSpeciesParameters& sP = this->speciesParams[popID];
+               for(uint i=0; i<6; i++) {
+                  if(isThisCellOnAFace[i] && facesToProcess[i] && !sP.facesToSkipVlasov[i]) {
+                     switch(sP.faceVlasovScheme[i]) {
+                        case vlasovscheme::NONE:
+                           break;
+                        case vlasovscheme::COPY:
+                           if(mpiGrid[cellID]->sysBoundaryLayer == 2) {
+                              const OutflowSpeciesParameters& sP = this->speciesParams[popID];
+                              vlasovBoundaryCopyFromTheClosestL1OutflowNbr(mpiGrid,cellID,true,popID,true); // first true means copy moments only, second true means V moments
+                              vlasovBoundaryCopyFromTheClosestL1OutflowNbr(mpiGrid,cellID,true,popID,false); // first true means copy moments only, second false means R moments
+                           }
+                           break;
+                        default:
+                           abort_mpi("ERROR: invalid Outflow Vlasov scheme", 1);
+                     } // switch
+                  } // if on face
+               } // faces
+            } // populations
+         } // cells
+      } // pragma omp parallel
+   }
+
    void Outflow::getFaces(bool* faces) {
       for(uint i=0; i<6; i++) {
          faces[i] = facesToProcess[i];

sysboundary/outflow.h

@@ -131,7 +131,8 @@ namespace SBC {
          const uint popID,
          const bool calculate_V_moments
       );
-
+      virtual void setupL2OutflowAtRestart(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid);
+
       virtual void getFaces(bool* faces);
       virtual std::string getName() const;
       virtual uint getIndex() const;

sysboundary/setmaxwellian.cpp

@@ -207,6 +207,20 @@ namespace SBC {
       B[2] = Bz;
 
       calculateCellMoments(&templateCell,true,false,true);
+      for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
+         Population &pop = templateCell.get_population(popID);
+         pop.RHO_V = pop.RHO;
+         pop.RHO_R = pop.RHO;
+         for(int i=0; i<3; i++) {
+            pop.V_R[i] = pop.V[i];
+            pop.V_V[i] = pop.V[i];
+         }
+         for(int i=0; i<6; i++) {
+            pop.P_R[i] = pop.P[i];
+            pop.P_V[i] = pop.P[i];
+         }
+      }
+
       templateCell.parameters[CellParams::RHOM_R] = templateCell.parameters[CellParams::RHOM];
       templateCell.parameters[CellParams::VX_R] = templateCell.parameters[CellParams::VX];
       templateCell.parameters[CellParams::VY_R] = templateCell.parameters[CellParams::VY];

sysboundary/sysboundary.cpp

@@ -649,6 +649,20 @@ void SysBoundary::updateState(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometr
    }
 }
 
+
+void SysBoundary::setupL2OutflowAtRestart(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid) {
+   list<SBC::SysBoundaryCondition*>::iterator it;
+   for (it = sysBoundaries.begin(); it != sysBoundaries.end(); it++) {
+      if (Parameters::isRestart           // When restarting
+         && !(*it)->doApplyUponRestart() // When reapplication is not requested
+         && (*it)->getIndex() == sysboundarytype::OUTFLOW
+      ) {
+         (*it)->setupL2OutflowAtRestart(mpiGrid);
+      }
+   }
+}
+
+
 /*!\brief Apply the Vlasov system boundary conditions to all system boundary cells at time t.
  *
  * Loops through all SysBoundaryConditions and calls the corresponding vlasovBoundaryCondition() function.
@@ -672,26 +686,25 @@ void SysBoundary::applySysBoundaryVlasovConditions(
 
 /*Transfer along boundaries*/
 // First the small stuff without overlapping in an extended neighbourhood:
-// TODO This now communicates in the wider neighbourhood for both layers, could be reduced to smaller neighbourhood for layer 1, larger neighbourhood for layer 2.
    SpatialCell::set_mpi_transfer_type(Transfer::CELL_PARAMETERS | Transfer::POP_METADATA | Transfer::CELL_SYSBOUNDARYFLAG, true);
    mpiGrid.update_copies_of_remote_neighbors(Neighborhoods::SYSBOUNDARIES_EXTENDED);
 
    // Loop over existing particle species
    for (uint popID = 0; popID < getObjectWrapper().particleSpecies.size(); ++popID) {
       SpatialCell::setCommunicatedSpecies(popID);
-      // update lists in larger neighborhood
-      updateRemoteVelocityBlockLists(mpiGrid, popID, Neighborhoods::SYSBOUNDARIES_EXTENDED);
+      // update lists in neighborhood
+      updateRemoteVelocityBlockLists(mpiGrid, popID, Neighborhoods::SYSBOUNDARIES);
 
       // Then the block data in the reduced neighbourhood:
       phiprof::Timer commTimer {"Start comm of cell and block data", {"MPI"}};
       SpatialCell::set_mpi_transfer_type(Transfer::VEL_BLOCK_DATA, true);
-      mpiGrid.start_remote_neighbor_copy_updates(Neighborhoods::SYSBOUNDARIES_EXTENDED);
+      mpiGrid.start_remote_neighbor_copy_updates(Neighborhoods::SYSBOUNDARIES);
       commTimer.stop();
 
       phiprof::Timer computeInnerTimer {"Compute process inner cells"};
       // Compute Vlasov boundary condition on system boundary/process inner cells
       vector<CellID> localCells;
-      getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_not_on_process_boundary(Neighborhoods::SYSBOUNDARIES_EXTENDED), localCells);
+      getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_not_on_process_boundary(Neighborhoods::SYSBOUNDARIES), localCells);
 
 #pragma omp parallel for
       for (uint i = 0; i < localCells.size(); i++) {
@@ -709,13 +722,13 @@ void SysBoundary::applySysBoundaryVlasovConditions(
       computeInnerTimer.stop();
 
       phiprof::Timer waitimer {"Wait for receives", {"MPI", "Wait"}};
-      mpiGrid.wait_remote_neighbor_copy_updates(Neighborhoods::SYSBOUNDARIES_EXTENDED);
+      mpiGrid.wait_remote_neighbor_copy_updates(Neighborhoods::SYSBOUNDARIES);
       waitimer.stop();
 
       // Compute vlasov boundary on system boundary/process boundary cells
       phiprof::Timer computeBoundaryTimer {"Compute process boundary cells"};
       vector<CellID> boundaryCells;
-      getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_on_process_boundary(Neighborhoods::SYSBOUNDARIES_EXTENDED),
+      getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_on_process_boundary(Neighborhoods::SYSBOUNDARIES),
                           boundaryCells);
 #pragma omp parallel for
       for (uint i = 0; i < boundaryCells.size(); i++) {

sysboundary/sysboundary.h

@@ -84,6 +84,8 @@ class SysBoundary {
                     FsGrid<std::array<Real, fsgrids::bgbfield::N_BGB>, FS_STENCIL_WIDTH>& BgBGrid,
                     creal t);
    void applySysBoundaryVlasovConditions(dccrg::Dccrg<spatial_cell::SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid, creal& t, const bool calculate_V_moments);
+   void setupL2OutflowAtRestart(dccrg::Dccrg<SpatialCell, dccrg::Cartesian_Geometry>& mpiGrid);
+
    unsigned int size() const;
    SBC::SysBoundaryCondition* getSysBoundary(cuint sysBoundaryType) const;
    bool isAnyDynamic() const;