diff --git a/sysboundary/ionosphere.h b/sysboundary/ionosphere.h index 73f8b6b2d..df15686b8 100644 --- a/sysboundary/ionosphere.h +++ b/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 diff --git a/sysboundary/outflow.cpp b/sysboundary/outflow.cpp index c98b3d00a..20fc75141 100644 --- a/sysboundary/outflow.cpp +++ b/sysboundary/outflow.cpp @@ -59,12 +59,12 @@ namespace SBC { Readparameters::addComposing(pop + "_outflow.reapplyFaceUponRestart", "List of faces on which outflow boundary conditions are to be reapplied upon restart ([xyz][+-])."); Readparameters::addComposing(pop + "_outflow.face", "List of faces on which outflow boundary conditions are to be applied ([xyz][+-])."); - Readparameters::add(pop + "_outflow.vlasovScheme_face_x+", "Scheme to use on the face x+ (Copy, Limit, None)", defStr); - Readparameters::add(pop + "_outflow.vlasovScheme_face_x-", "Scheme to use on the face x- (Copy, Limit, None)", defStr); - Readparameters::add(pop + "_outflow.vlasovScheme_face_y+", "Scheme to use on the face y+ (Copy, Limit, None)", defStr); - Readparameters::add(pop + "_outflow.vlasovScheme_face_y-", "Scheme to use on the face y- (Copy, Limit, None)", defStr); - Readparameters::add(pop + "_outflow.vlasovScheme_face_z+", "Scheme to use on the face z+ (Copy, Limit, None)", defStr); - Readparameters::add(pop + "_outflow.vlasovScheme_face_z-", "Scheme to use on the face z- (Copy, Limit, None)", defStr); + Readparameters::add(pop + "_outflow.vlasovScheme_face_x+", "Scheme to use on the face x+ (Copy, None)", defStr); + Readparameters::add(pop + "_outflow.vlasovScheme_face_x-", "Scheme to use on the face x- (Copy, None)", defStr); + Readparameters::add(pop + "_outflow.vlasovScheme_face_y+", "Scheme to use on the face y+ (Copy, None)", defStr); + Readparameters::add(pop + "_outflow.vlasovScheme_face_y-", "Scheme to use on the face y- (Copy, None)", defStr); + Readparameters::add(pop + "_outflow.vlasovScheme_face_z+", "Scheme to use on the face z+ (Copy, None)", defStr); + Readparameters::add(pop + "_outflow.vlasovScheme_face_z-", "Scheme to use on the face z- (Copy, None)", defStr); Readparameters::add(pop + "_outflow.quench", "Factor by which to quench the inflowing parts of the velocity distribution function.", 1.0); } @@ -118,8 +118,6 @@ namespace SBC { sP.faceVlasovScheme[j] = vlasovscheme::NONE; } else if (vlasovSysBoundarySchemeName[j] == "Copy") { sP.faceVlasovScheme[j] = vlasovscheme::COPY; - } else if(vlasovSysBoundarySchemeName[j] == "Limit") { - sP.faceVlasovScheme[j] = vlasovscheme::LIMIT; } else { abort_mpi("ERROR: " + vlasovSysBoundarySchemeName[j] + " is an invalid Outflow Vlasov scheme!"); } @@ -271,7 +269,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]; @@ -415,10 +413,11 @@ namespace SBC { case vlasovscheme::NONE: break; case vlasovscheme::COPY: - vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,false,popID,calculate_V_moments); - break; - case vlasovscheme::LIMIT: - vlasovBoundaryCopyFromTheClosestNbrAndLimit(mpiGrid,cellID,popID); + 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); diff --git a/sysboundary/outflow.h b/sysboundary/outflow.h index a947d748b..fdb0015fa 100644 --- a/sysboundary/outflow.h +++ b/sysboundary/outflow.h @@ -153,7 +153,6 @@ namespace SBC { enum vlasovscheme { NONE, COPY, - LIMIT, N_SCHEMES }; }; // class Outflow diff --git a/sysboundary/sysboundary.cpp b/sysboundary/sysboundary.cpp index 6b355f3dd..275b9d400 100644 --- a/sysboundary/sysboundary.cpp +++ b/sysboundary/sysboundary.cpp @@ -663,26 +663,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(SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID); // 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, SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID); + // update lists in neighborhood + updateRemoteVelocityBlockLists(mpiGrid, popID, SYSBOUNDARIES_NEIGHBORHOOD_ID); // 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(SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID); + mpiGrid.start_remote_neighbor_copy_updates(SYSBOUNDARIES_NEIGHBORHOOD_ID); commTimer.stop(); phiprof::Timer computeInnerTimer {"Compute process inner cells"}; // Compute Vlasov boundary condition on system boundary/process inner cells vector localCells; - getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_not_on_process_boundary(SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID), localCells); + getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_not_on_process_boundary(SYSBOUNDARIES_NEIGHBORHOOD_ID), localCells); #pragma omp parallel for for (uint i = 0; i < localCells.size(); i++) { @@ -700,13 +699,13 @@ void SysBoundary::applySysBoundaryVlasovConditions( computeInnerTimer.stop(); phiprof::Timer waitimer {"Wait for receives", {"MPI", "Wait"}}; - mpiGrid.wait_remote_neighbor_copy_updates(SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID); + mpiGrid.wait_remote_neighbor_copy_updates(SYSBOUNDARIES_NEIGHBORHOOD_ID); waitimer.stop(); // Compute vlasov boundary on system boundary/process boundary cells phiprof::Timer computeBoundaryTimer {"Compute process boundary cells"}; vector boundaryCells; - getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_on_process_boundary(SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID), + getBoundaryCellList(mpiGrid, mpiGrid.get_local_cells_on_process_boundary(SYSBOUNDARIES_NEIGHBORHOOD_ID), boundaryCells); #pragma omp parallel for for (uint i = 0; i < boundaryCells.size(); i++) { diff --git a/sysboundary/sysboundarycondition.cpp b/sysboundary/sysboundarycondition.cpp index e39d9927a..b3efa9dbc 100644 --- a/sysboundary/sysboundarycondition.cpp +++ b/sysboundary/sysboundarycondition.cpp @@ -341,78 +341,6 @@ namespace SBC { averageCellData(mpiGrid, closeCells, mpiGrid[cellID], popID, fluffiness); } - /*! Function used to copy the distribution from (one of) the closest sysboundarytype::NOT_SYSBOUNDARY cell but limiting to values no higher than where it can flow into. Moments are recomputed. - * \param mpiGrid Grid - * \param cellID The cell's ID. - */ - void SysBoundaryCondition::vlasovBoundaryCopyFromTheClosestNbrAndLimit( - dccrg::Dccrg& mpiGrid, - const CellID& cellID, - const uint popID - ) { - const CellID closestCell = getTheClosestNonsysboundaryCell(cellID); - SpatialCell * from = mpiGrid[closestCell]; - SpatialCell * to = mpiGrid[cellID]; - - if(closestCell == INVALID_CELLID) { - abort_mpi("No closest cell found!", 1); - } - - const array& flowtoCells = getFlowtoCells(cellID); - //Do not allow block adjustment, the block structure when calling vlasovBoundaryCondition should be static - //just copy data to existing blocks, no modification of to blocks allowed - for (vmesh::LocalID blockLID=0; blockLIDget_number_of_velocity_blocks(popID); ++blockLID) { - const vmesh::GlobalID blockGID = to->get_velocity_block_global_id(blockLID,popID); -// const Realf* fromBlock_data = from->get_data(from->get_velocity_block_local_id(blockGID) ); - Realf* toBlock_data = to->get_data(blockLID,popID); - if (from->get_velocity_block_local_id(blockGID,popID) == from->invalid_local_id()) { - for (unsigned int i = 0; i < VELOCITY_BLOCK_LENGTH; i++) { - toBlock_data[i] = 0.0; //block did not exist in from cell, fill with zeros. - } - } else { - const Real* blockParameters = to->get_block_parameters(blockLID, popID); - // check where cells are - creal vxBlock = blockParameters[BlockParams::VXCRD]; - creal vyBlock = blockParameters[BlockParams::VYCRD]; - creal vzBlock = blockParameters[BlockParams::VZCRD]; - creal dvxCell = blockParameters[BlockParams::DVX]; - creal dvyCell = blockParameters[BlockParams::DVY]; - creal dvzCell = blockParameters[BlockParams::DVZ]; - - array flowtoCellsBlockCache = getFlowtoCellsBlock(flowtoCells, blockGID, popID); - - for (uint kc=0; kcget_velocity_cell(indices); - - creal vxCellCenter = vxBlock + (ic+convert(0.5))*dvxCell; - creal vyCellCenter = vyBlock + (jc+convert(0.5))*dvyCell; - creal vzCellCenter = vzBlock + (kc+convert(0.5))*dvzCell; - const int vxCellSign = vxCellCenter < 0 ? -1 : 1; - const int vyCellSign = vyCellCenter < 0 ? -1 : 1; - const int vzCellSign = vzCellCenter < 0 ? -1 : 1; - Realf value = from->get_value(blockGID, cell, popID); - //loop over spatial cells in quadrant of influence - for(int dvx = 0 ; dvx <= 1; dvx++) { - for(int dvy = 0 ; dvy <= 1; dvy++) { - for(int dvz = 0 ; dvz <= 1; dvz++) { - const int flowToId = nbrID(dvx * vxCellSign, dvy * vyCellSign, dvz * vzCellSign); - if(flowtoCells.at(flowToId)){ - value = min(value, flowtoCellsBlockCache.at(flowToId)[cell]); - } - } - } - } - to->set_value(blockGID, cell, value, popID); - } - } - } - } - } - } - /*! Function used to copy the distribution and moments from one cell to another. * \param from Pointer to parent cell to copy from. * \param to Pointer to destination cell. @@ -446,9 +374,7 @@ namespace SBC { } } - if(!copyMomentsOnly) { // Do this only if copyMomentsOnly is false. - to->set_population(from->get_population(popID), popID); - } else { + if(copyMomentsOnly) { if (calculate_V_moments) { to->get_population(popID).RHO_V = from->get_population(popID).RHO_V; } else { @@ -464,6 +390,8 @@ namespace SBC { to->get_population(popID).P_R[i] = from->get_population(popID).P_R[i]; } } + } else { + to->set_population(from->get_population(popID), popID); } } diff --git a/sysboundary/sysboundarycondition.h b/sysboundary/sysboundarycondition.h index b06592e7d..6f6bfc1b6 100644 --- a/sysboundary/sysboundarycondition.h +++ b/sysboundary/sysboundarycondition.h @@ -227,11 +227,6 @@ namespace SBC { const uint popID, const bool calculate_V_moments ); - void vlasovBoundaryCopyFromTheClosestNbrAndLimit( - dccrg::Dccrg& mpiGrid, - const CellID& cellID, - const uint popID - ); void vlasovBoundaryCopyFromAllClosestNbrs( dccrg::Dccrg& mpiGrid, const CellID& cellID, diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 7e032674a..bae8ed069 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -165,9 +165,11 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, * @param cells Vector containing the spatial cells to be calculated. * @param computeSecond If true, second velocity moments are calculated.*/ void calculateMoments_R( - dccrg::Dccrg& mpiGrid, - const std::vector& cells, - const bool& computeSecond) { + dccrg::Dccrg& mpiGrid, + const std::vector& cells, + const bool& computeSecond, + const bool initialCompute +) { phiprof::Timer momentsTimer {"compute-moments-n"}; @@ -179,6 +181,9 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code + continue; + } // Clear old moments to zero value if (popID == 0) { @@ -255,6 +260,9 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code + continue; + } cell->parameters[CellParams::VX_R] = divideIfNonZero(cell->parameters[CellParams::VX_R], cell->parameters[CellParams::RHOM_R]); cell->parameters[CellParams::VY_R] = divideIfNonZero(cell->parameters[CellParams::VY_R], cell->parameters[CellParams::RHOM_R]); cell->parameters[CellParams::VZ_R] = divideIfNonZero(cell->parameters[CellParams::VZ_R], cell->parameters[CellParams::RHOM_R]); @@ -273,6 +281,9 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code + continue; + } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); if (blockContainer.size() == 0) { @@ -321,9 +332,11 @@ void calculateMoments_R( * @param cells Vector containing the spatial cells to be calculated. * @param computeSecond If true, second velocity moments are calculated.*/ void calculateMoments_V( - dccrg::Dccrg& mpiGrid, - const std::vector& cells, - const bool& computeSecond) { + dccrg::Dccrg& mpiGrid, + const std::vector& cells, + const bool& computeSecond, + const bool initialCompute +) { phiprof::Timer momentsTimer {"Compute _V moments"}; @@ -336,6 +349,9 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code + continue; + } // Clear old moments to zero value if (popID == 0) { @@ -398,6 +414,10 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code + continue; + } + cell->parameters[CellParams::VX_V] = divideIfNonZero(cell->parameters[CellParams::VX_V], cell->parameters[CellParams::RHOM_V]); cell->parameters[CellParams::VY_V] = divideIfNonZero(cell->parameters[CellParams::VY_V], cell->parameters[CellParams::RHOM_V]); cell->parameters[CellParams::VZ_V] = divideIfNonZero(cell->parameters[CellParams::VZ_V], cell->parameters[CellParams::RHOM_V]); @@ -416,6 +436,9 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code + continue; + } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); if (blockContainer.size() == 0) { diff --git a/vlasovsolver/cpu_moments.h b/vlasovsolver/cpu_moments.h index 962920c41..630e3bb15 100644 --- a/vlasovsolver/cpu_moments.h +++ b/vlasovsolver/cpu_moments.h @@ -45,13 +45,19 @@ void blockVelocitySecondMoments(const Realf* avgs,const Real* blockParams, const REAL v[3], std::vector& array); -void calculateMoments_R(dccrg::Dccrg& mpiGrid, - const std::vector& cells, - const bool& computeSecond); - -void calculateMoments_V(dccrg::Dccrg& mpiGrid, - const std::vector& cells, - const bool& computeSecond); +void calculateMoments_R( + dccrg::Dccrg& mpiGrid, + const std::vector& cells, + const bool& computeSecond, + const bool initialCompute=false +); + +void calculateMoments_V( + dccrg::Dccrg& mpiGrid, + const std::vector& cells, + const bool& computeSecond, + const bool initialCompute=false +); diff --git a/vlasovsolver/cpu_trans_pencils.cpp b/vlasovsolver/cpu_trans_pencils.cpp index 38cfc40ee..589837f4d 100644 --- a/vlasovsolver/cpu_trans_pencils.cpp +++ b/vlasovsolver/cpu_trans_pencils.cpp @@ -295,7 +295,9 @@ void computeSpatialSourceCellsForPencil(const dccrg::DccrgsysBoundaryFlag != sysboundarytype::DO_NOT_COMPUTE) { + if (mpiGrid[ids[i]]->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY || + (mpiGrid[ids[i]]->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && mpiGrid[ids[i]]->sysBoundaryFlag != sysboundarytype::DO_NOT_COMPUTE && mpiGrid[ids[i]]->sysBoundaryLayer == 1) + ) { isGood = true; } } @@ -313,7 +315,9 @@ void computeSpatialSourceCellsForPencil(const dccrg::DccrgsysBoundaryFlag != sysboundarytype::DO_NOT_COMPUTE) { + if (mpiGrid[ids[i]]->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY || + (mpiGrid[ids[i]]->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && mpiGrid[ids[i]]->sysBoundaryFlag != sysboundarytype::DO_NOT_COMPUTE && mpiGrid[ids[i]]->sysBoundaryLayer == 1) + ) { isGood = true; } } diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 62e290a9c..32bd85dda 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -257,7 +257,7 @@ void calculateSpatialTranslation( // If dt=0 we are either initializing or distribution functions are not translated. // In both cases go to the end of this function and calculate the moments. if (dt == 0.0) { - calculateMoments_R(mpiGrid,localCells,true); + calculateMoments_R(mpiGrid,localCells,true,true); return; } @@ -493,7 +493,7 @@ void calculateAcceleration(dccrg::Dccrg& } // Recalculate "_V" velocity moments - calculateMoments_V(mpiGrid,cells,true); + calculateMoments_V(mpiGrid,cells,true,(dt==0)); // Set CellParams::MAXVDT to be the minimum dt of all per-species values #pragma omp parallel for