Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SysBoundary communicate VDFs only in smaller neighborhood #1049

Open
wants to merge 12 commits into
base: dev
Choose a base branch
from
Open

SysBoundary communicate VDFs only in smaller neighborhood #1049

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
7e65124
SysBoundary communicate VDFs only in smaller neighborhood
ykempf Oct 31, 2024
beb7333
Merge remote-tracking branch 'fmihpc/dev' into less_L2_boundary_vdf_c…
ykempf Nov 4, 2024
5b8b97c
No copy of VDFs for L2 outflow cells.
ykempf Nov 4, 2024
550424b
Remove LIMIT outflow scheme to avoid L2 VDF requirements.
ykempf Nov 4, 2024
d086209
Updated comments in sysboundary about the change in communication width.
ykempf Nov 4, 2024
b765acc
Removed fluffiness parameter from ionosphere.h.
ykempf Nov 4, 2024
0f13e82
Try computing _R and _V moments at t=0 for everyone.
ykempf Nov 1, 2024
826f681
Try to get moments correct at restart too.
ykempf Nov 5, 2024
fdfc499
Change order of if / else in copyCellData
ykempf Nov 12, 2024
ccd369a
Fix initialization of _R and _V moments in L2 boundary cells.
ykempf Nov 12, 2024
83e40e7
Merge remote-tracking branch 'fmihpc/dev' into less_L2_boundary_vdf_c…
ykempf Nov 12, 2024
c8bad25
Flag non-L1 boundary cells as not good for translation pencils.
ykempf Nov 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion sysboundary/ionosphere.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
25 changes: 12 additions & 13 deletions sysboundary/outflow.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand Down Expand Up @@ -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!");
}
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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);
Expand Down
1 change: 0 additions & 1 deletion sysboundary/outflow.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,6 @@ namespace SBC {
enum vlasovscheme {
NONE,
COPY,
LIMIT,
N_SCHEMES
};
}; // class Outflow
Expand Down
13 changes: 6 additions & 7 deletions sysboundary/sysboundary.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<CellID> 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++) {
Expand All @@ -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<CellID> 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++) {
Expand Down
78 changes: 3 additions & 75 deletions sysboundary/sysboundarycondition.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<SpatialCell,dccrg::Cartesian_Geometry>& 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<SpatialCell*,27>& 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; blockLID<to->get_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<Realf*,27> flowtoCellsBlockCache = getFlowtoCellsBlock(flowtoCells, blockGID, popID);

for (uint kc=0; kc<WID; ++kc) {
for (uint jc=0; jc<WID; ++jc) {
for (uint ic=0; ic<WID; ++ic) {
velocity_cell_indices_t indices = {ic, jc, kc};
const uint cell = from->get_velocity_cell(indices);

creal vxCellCenter = vxBlock + (ic+convert<Real>(0.5))*dvxCell;
creal vyCellCenter = vyBlock + (jc+convert<Real>(0.5))*dvyCell;
creal vzCellCenter = vzBlock + (kc+convert<Real>(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.
Expand Down Expand Up @@ -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 {
Expand All @@ -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);
}
}

Expand Down
5 changes: 0 additions & 5 deletions sysboundary/sysboundarycondition.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,11 +227,6 @@ namespace SBC {
const uint popID,
const bool calculate_V_moments
);
void vlasovBoundaryCopyFromTheClosestNbrAndLimit(
dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const CellID& cellID,
const uint popID
);
void vlasovBoundaryCopyFromAllClosestNbrs(
dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const CellID& cellID,
Expand Down
35 changes: 29 additions & 6 deletions vlasovsolver/cpu_moments.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond) {
dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond,
const bool initialCompute
) {

phiprof::Timer momentsTimer {"compute-moments-n"};

Expand All @@ -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) {
Expand Down Expand Up @@ -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]);
Expand All @@ -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<vmesh::LocalID>& blockContainer = cell->get_velocity_blocks(popID);
if (blockContainer.size() == 0) {
Expand Down Expand Up @@ -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<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond) {
dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond,
const bool initialCompute
) {

phiprof::Timer momentsTimer {"Compute _V moments"};

Expand All @@ -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) {
Expand Down Expand Up @@ -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]);
Expand All @@ -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<vmesh::LocalID>& blockContainer = cell->get_velocity_blocks(popID);
if (blockContainer.size() == 0) {
Expand Down
20 changes: 13 additions & 7 deletions vlasovsolver/cpu_moments.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,19 @@ void blockVelocitySecondMoments(const Realf* avgs,const Real* blockParams,
const REAL v[3],
std::vector<Real>& array);

void calculateMoments_R(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond);

void calculateMoments_V(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond);
void calculateMoments_R(
dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond,
const bool initialCompute=false
);

void calculateMoments_V(
dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
const std::vector<CellID>& cells,
const bool& computeSecond,
const bool initialCompute=false
);



Expand Down
8 changes: 6 additions & 2 deletions vlasovsolver/cpu_trans_pencils.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,9 @@ void computeSpatialSourceCellsForPencil(const dccrg::Dccrg<SpatialCell,dccrg::Ca
bool isGood = false;
if (ids[i]!=0) {
if (mpiGrid[ids[i]] != NULL) {
if (mpiGrid[ids[i]]->sysBoundaryFlag != 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;
}
}
Expand All @@ -313,7 +315,9 @@ void computeSpatialSourceCellsForPencil(const dccrg::Dccrg<SpatialCell,dccrg::Ca
bool isGood = false;
if (ids[i]!=0) {
if (mpiGrid[ids[i]] != NULL) {
if (mpiGrid[ids[i]]->sysBoundaryFlag != 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;
}
}
Expand Down
4 changes: 2 additions & 2 deletions vlasovsolver/vlasovmover.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand Down Expand Up @@ -493,7 +493,7 @@ void calculateAcceleration(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>&
}

// 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
Expand Down