From 7e6512478870b4513bf27ec1e4b532f867524253 Mon Sep 17 00:00:00 2001 From: ykempf Date: Thu, 31 Oct 2024 16:33:15 +0200 Subject: [PATCH 01/10] SysBoundary communicate VDFs only in smaller neighborhood --- sysboundary/sysboundary.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sysboundary/sysboundary.cpp b/sysboundary/sysboundary.cpp index 008ba8417..db3032335 100644 --- a/sysboundary/sysboundary.cpp +++ b/sysboundary/sysboundary.cpp @@ -677,18 +677,18 @@ void SysBoundary::applySysBoundaryVlasovConditions( for (uint popID = 0; popID < getObjectWrapper().particleSpecies.size(); ++popID) { SpatialCell::setCommunicatedSpecies(popID); // update lists in larger neighborhood - updateRemoteVelocityBlockLists(mpiGrid, popID, SYSBOUNDARIES_EXTENDED_NEIGHBORHOOD_ID); + 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++) { @@ -706,13 +706,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++) { From 5b8b97c989311aa122c24cbeba6849a27e7f890a Mon Sep 17 00:00:00 2001 From: ykempf Date: Mon, 4 Nov 2024 18:06:12 +0200 Subject: [PATCH 02/10] No copy of VDFs for L2 outflow cells. --- sysboundary/outflow.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/sysboundary/outflow.cpp b/sysboundary/outflow.cpp index c98b3d00a..97ce3e2cf 100644 --- a/sysboundary/outflow.cpp +++ b/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]; @@ -415,10 +415,14 @@ 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; case vlasovscheme::LIMIT: - vlasovBoundaryCopyFromTheClosestNbrAndLimit(mpiGrid,cellID,popID); + vlasovBoundaryCopyFromTheClosestNbrAndLimit(mpiGrid,cellID,popID); // may do VDF but easier this way break; default: abort_mpi("ERROR: invalid Outflow Vlasov scheme", 1); From 550424bc15a68a31c11ffcdd7110440be74da4ce Mon Sep 17 00:00:00 2001 From: ykempf Date: Mon, 4 Nov 2024 18:24:54 +0200 Subject: [PATCH 03/10] Remove LIMIT outflow scheme to avoid L2 VDF requirements. --- sysboundary/outflow.cpp | 17 +++---- sysboundary/outflow.h | 1 - sysboundary/sysboundarycondition.cpp | 74 +--------------------------- sysboundary/sysboundarycondition.h | 5 -- 4 files changed, 7 insertions(+), 90 deletions(-) diff --git a/sysboundary/outflow.cpp b/sysboundary/outflow.cpp index 97ce3e2cf..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!"); } @@ -421,9 +419,6 @@ namespace SBC { vlasovBoundaryCopyFromTheClosestNbr(mpiGrid,cellID,true,popID,calculate_V_moments); // no VDF copy } break; - case vlasovscheme::LIMIT: - vlasovBoundaryCopyFromTheClosestNbrAndLimit(mpiGrid,cellID,popID); // may do VDF but easier this way - 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/sysboundarycondition.cpp b/sysboundary/sysboundarycondition.cpp index e39d9927a..ff4eb0131 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,7 +374,7 @@ namespace SBC { } } - if(!copyMomentsOnly) { // Do this only if copyMomentsOnly is false. + if(!copyMomentsOnly) { // Do this only if copyMomentsOnly is false. to->set_population(from->get_population(popID), popID); } else { if (calculate_V_moments) { 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, From d0862093048ce40dd44c6dd4816314f761d05fc1 Mon Sep 17 00:00:00 2001 From: ykempf Date: Mon, 4 Nov 2024 18:31:24 +0200 Subject: [PATCH 04/10] Updated comments in sysboundary about the change in communication width. --- sysboundary/sysboundary.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sysboundary/sysboundary.cpp b/sysboundary/sysboundary.cpp index 9839fc7bd..275b9d400 100644 --- a/sysboundary/sysboundary.cpp +++ b/sysboundary/sysboundary.cpp @@ -663,14 +663,13 @@ 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 + // update lists in neighborhood updateRemoteVelocityBlockLists(mpiGrid, popID, SYSBOUNDARIES_NEIGHBORHOOD_ID); // Then the block data in the reduced neighbourhood: From b765acc7d410afc2a6fa2e6010caf6fdd74b55e9 Mon Sep 17 00:00:00 2001 From: ykempf Date: Mon, 4 Nov 2024 18:45:17 +0200 Subject: [PATCH 05/10] Removed fluffiness parameter from ionosphere.h. --- sysboundary/ionosphere.h | 1 - 1 file changed, 1 deletion(-) 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 From 0f13e820810a6f105b61d1c89fd601edc9d6902e Mon Sep 17 00:00:00 2001 From: ykempf Date: Fri, 1 Nov 2024 15:12:34 +0200 Subject: [PATCH 06/10] Try computing _R and _V moments at t=0 for everyone. --- vlasovsolver/cpu_moments.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 995dc12f5..61b4c6662 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -176,6 +176,9 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + continue; + } // Clear old moments to zero value if (popID == 0) { @@ -249,6 +252,9 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // 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]); @@ -267,6 +273,9 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + continue; + } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); if (blockContainer.size() == 0) { @@ -330,6 +339,9 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + continue; + } // Clear old moments to zero value if (popID == 0) { @@ -389,6 +401,10 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // 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]); @@ -407,6 +423,9 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + continue; + } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); if (blockContainer.size() == 0) { From 826f681e150bde0fc19c8564884f62948d2d2f9c Mon Sep 17 00:00:00 2001 From: ykempf Date: Tue, 5 Nov 2024 13:03:40 +0200 Subject: [PATCH 07/10] Try to get moments correct at restart too. --- vlasovsolver/cpu_moments.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 61b4c6662..d0d46fc0c 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -176,7 +176,7 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code continue; } @@ -252,7 +252,7 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // 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]); @@ -273,7 +273,7 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code continue; } @@ -339,7 +339,7 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code continue; } @@ -401,7 +401,7 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code continue; } @@ -423,7 +423,7 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != 0) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code continue; } From fdfc4999b1186cfb3b178856a1540588ebfddc92 Mon Sep 17 00:00:00 2001 From: ykempf Date: Tue, 12 Nov 2024 12:17:47 +0200 Subject: [PATCH 08/10] Change order of if / else in copyCellData --- sysboundary/sysboundarycondition.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sysboundary/sysboundarycondition.cpp b/sysboundary/sysboundarycondition.cpp index ff4eb0131..b3efa9dbc 100644 --- a/sysboundary/sysboundarycondition.cpp +++ b/sysboundary/sysboundarycondition.cpp @@ -374,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 { @@ -392,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); } } From ccd369ae4c5152fb3b64005b8bf99a1712142c26 Mon Sep 17 00:00:00 2001 From: ykempf Date: Tue, 12 Nov 2024 12:18:55 +0200 Subject: [PATCH 09/10] Fix initialization of _R and _V moments in L2 boundary cells. --- vlasovsolver/cpu_moments.cpp | 28 ++++++++++++++++------------ vlasovsolver/cpu_moments.h | 20 +++++++++++++------- vlasovsolver/vlasovmover.cpp | 4 ++-- 3 files changed, 31 insertions(+), 21 deletions(-) diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index d0d46fc0c..ef6792c1c 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -162,9 +162,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"}; @@ -176,7 +178,7 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code continue; } @@ -252,7 +254,7 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code + 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]); @@ -273,7 +275,7 @@ void calculateMoments_R( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code continue; } @@ -324,9 +326,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"}; @@ -339,7 +343,7 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code continue; } @@ -401,7 +405,7 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code continue; } @@ -423,7 +427,7 @@ void calculateMoments_V( if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) { continue; } - if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && P::tstep != P::tstep_min) { // these should have been handled by the boundary code + if (cell->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY && cell->sysBoundaryLayer != 1 && !initialCompute) { // these should have been handled by the boundary code continue; } diff --git a/vlasovsolver/cpu_moments.h b/vlasovsolver/cpu_moments.h index ce98931cb..4d8c06931 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], REAL* 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/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 0332df94f..cbf421a83 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 From c8bad25695cbf84d56bf19d944d8341ad34bd5cb Mon Sep 17 00:00:00 2001 From: ykempf Date: Thu, 28 Nov 2024 11:38:21 +0200 Subject: [PATCH 10/10] Flag non-L1 boundary cells as not good for translation pencils. --- vlasovsolver/cpu_trans_pencils.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) 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; } }