diff --git a/common.h b/common.h index ef76de4e3..253f3df4c 100644 --- a/common.h +++ b/common.h @@ -165,15 +165,27 @@ namespace CellParams { P_11, /*!< Pressure P_xx component, computed by Vlasov propagator. */ P_22, /*!< Pressure P_yy component, computed by Vlasov propagator. */ P_33, /*!< Pressure P_zz component, computed by Vlasov propagator. */ + P_23, /*!< Pressure P_yz component, computed by Vlasov propagator. */ + P_13, /*!< Pressure P_xz component, computed by Vlasov propagator. */ + P_12, /*!< Pressure P_xy component, computed by Vlasov propagator. */ P_11_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_11_R and P_11_V. */ P_22_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_22_R and P_22_V. */ P_33_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_33_R and P_33_V. */ + P_23_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_23_R and P_23_V. */ + P_13_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_13_R and P_13_V. */ + P_12_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_12_R and P_12_V. */ P_11_R, /*!< P_xx component after propagation in ordinary space */ P_22_R, /*!< P_yy component after propagation in ordinary space */ P_33_R, /*!< P_zz component after propagation in ordinary space */ + P_23_R, /*!< P_yz component after propagation in ordinary space */ + P_13_R, /*!< P_xz component after propagation in ordinary space */ + P_12_R, /*!< P_xy component after propagation in ordinary space */ P_11_V, /*!< P_xx component after propagation in velocity space */ P_22_V, /*!< P_yy component after propagation in velocity space */ P_33_V, /*!< P_zz component after propagation in velocity space */ + P_23_V, /*!< P_yz component after propagation in velocity space */ + P_12_V, /*!< P_xy component after propagation in velocity space */ + P_13_V, /*!< P_xz component after propagation in velocity space */ EXVOL, /*!< Volume electric field averaged over spatial cell, x-component.*/ EYVOL, /*!< Volume electric field averaged over spatial cell, y-component.*/ EZVOL, /*!< Volume electric field averaged over spatial cell, z-component.*/ @@ -209,6 +221,8 @@ namespace CellParams { AMR_DB, AMR_ALPHA1, AMR_ALPHA2, + P_ANISOTROPY, + AMR_VORTICITY, RECENTLY_REFINED, BULKV_FORCING_X, /*! Externally forced drift velocity (ex. from the ionosphere) */ BULKV_FORCING_Y, /*! Externally forced drift velocity (ex. from the ionosphere) */ @@ -235,6 +249,22 @@ namespace bvolderivatives { }; } +namespace vderivatives { + // Essentially a copy from dmoments + enum { + dVxdx, /*!< Derivative of volume-averaged Vx to x-direction. */ + dVxdy, /*!< Derivative of volume-averaged Vx to y-direction. */ + dVxdz, /*!< Derivative of volume-averaged Vx to z-direction. */ + dVydx, /*!< Derivative of volume-averaged Vy to x-direction. */ + dVydy, /*!< Derivative of volume-averaged Vy to y-direction. */ + dVydz, /*!< Derivative of volume-averaged Vy to z-direction. */ + dVzdx, /*!< Derivative of volume-averaged Vz to x-direction. */ + dVzdy, /*!< Derivative of volume-averaged Vz to y-direction. */ + dVzdz, /*!< Derivative of volume-averaged Vz to z-direction. */ + N_V_DERIVATIVES + }; +} + // FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, // FsGrid< std::array, FS_STENCIL_WIDTH> & perBDt2Grid, // FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp index 309915ae9..dda5607d0 100644 --- a/datareduction/datareducer.cpp +++ b/datareduction/datareducer.cpp @@ -2962,7 +2962,21 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti } if(P::systemWriteAllDROs || lowercase == "vg_amr_alpha2") { outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_amr_alpha2",CellParams::AMR_ALPHA2,1)); - outputReducer->addMetadata(outputReducer->size()-1,"","","$\\alpha_2",""); + outputReducer->addMetadata(outputReducer->size()-1,"","","$\\alpha_2$",""); + if(!P::systemWriteAllDROs) { + continue; + } + } + if(P::systemWriteAllDROs || lowercase == "vg_pressure_anisotropy") { + outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_pressure_anisotropy",CellParams::P_ANISOTROPY,1)); + outputReducer->addMetadata(outputReducer->size()-1,"","","$P_\\perp / P_\\parallel$",""); + if(!P::systemWriteAllDROs) { + continue; + } + } + if(P::systemWriteAllDROs || lowercase == "vg_amr_vorticity") { + outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_amr_vorticity",CellParams::AMR_VORTICITY,1)); + outputReducer->addMetadata(outputReducer->size()-1,"","","Vorticity",""); if(!P::systemWriteAllDROs) { continue; } diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index 5d3bc3ca6..8112d08ac 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -25,6 +25,7 @@ #include "fs_common.h" #include "derivatives.hpp" #include "fs_limiters.h" +#include /*! \brief Low-level spatial derivatives calculation. * @@ -758,15 +759,22 @@ void calculateScaledDeltas( Real myRho {cell->parameters[CellParams::RHOM]}; Real myU {calculateU(cell)}; + Real myV {std::sqrt(std::pow(cell->parameters[CellParams::VX], 2) + std::pow(cell->parameters[CellParams::VY], 2) + std::pow(cell->parameters[CellParams::VZ], 2))}; + Real maxV {myV}; std::array myP = getMomentumDensity(cell); std::array myB = getBVol(cell); for (SpatialCell* neighbor : neighbors) { Real otherRho = neighbor->parameters[CellParams::RHOM]; Real otherU = calculateU(neighbor); + Real otherV {std::sqrt(std::pow(neighbor->parameters[CellParams::VX], 2) + std::pow(neighbor->parameters[CellParams::VY], 2) + std::pow(neighbor->parameters[CellParams::VZ], 2))}; std::array otherP = getMomentumDensity(neighbor); std::array otherB = getBVol(neighbor); Real deltaBsq = pow(myB[0] - otherB[0], 2) + pow(myB[1] - otherB[1], 2) + pow(myB[2] - otherB[2], 2); + if (myV < EPS) { + maxV = std::max(maxV,otherV); + } + Real maxRho = std::max(myRho, otherRho); if (maxRho > EPS) { dRho = std::max(fabs(myRho - otherRho) / maxRho, dRho); @@ -820,6 +828,38 @@ void calculateScaledDeltas( Bperp = std::sqrt(Bperp); } + // Now, rotation matrix to get parallel and perpendicular pressure + //Eigen::Quaterniond q {Quaterniond::FromTwoVectors(Eigen::vector3d{0, 0, 1}, Eigen::vector3d{myB[0], myB[1], myB[2]})}; + //Eigen::Matrix3d rot = q.toRotationMatrix(); + Eigen::Matrix3d rot = Eigen::Quaterniond::FromTwoVectors(Eigen::Vector3d{myB[0], myB[1], myB[2]}, Eigen::Vector3d{0, 0, 1}).normalized().toRotationMatrix(); + Eigen::Matrix3d P { + {cell->parameters[CellParams::P_11], cell->parameters[CellParams::P_12], cell->parameters[CellParams::P_13]}, + {cell->parameters[CellParams::P_12], cell->parameters[CellParams::P_22], cell->parameters[CellParams::P_23]}, + {cell->parameters[CellParams::P_13], cell->parameters[CellParams::P_23], cell->parameters[CellParams::P_33]}, + }; + + Eigen::Matrix3d transposerot = rot.transpose(); + Eigen::Matrix3d Pprime = rot * P * transposerot; + + Real Panisotropy {0.0}; + if (Pprime(2, 2) > EPS) { + Panisotropy = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); + } + + // Vorticity + Real dVxdy {cell->derivativesV[vderivatives::dVxdy]}; + Real dVxdz {cell->derivativesV[vderivatives::dVxdz]}; + Real dVydx {cell->derivativesV[vderivatives::dVydx]}; + Real dVydz {cell->derivativesV[vderivatives::dVydz]}; + Real dVzdx {cell->derivativesV[vderivatives::dVzdx]}; + Real dVzdy {cell->derivativesV[vderivatives::dVzdy]}; + Real vorticity {std::sqrt(std::pow(dVxdy - dVydz, 2) + std::pow(dVxdz - dVzdx, 2 ) + std::pow(dVydx - dVxdy, 2))}; + //Real vA {std::sqrt(Bsq / (physicalconstants::MU_0 * myRho))}; + Real amr_vorticity {0.0}; + if (maxV > EPS) { + amr_vorticity = vorticity * cell->parameters[CellParams::DX] / maxV; + } + cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; cell->parameters[CellParams::AMR_DPSQ] = dPsq; @@ -827,6 +867,9 @@ void calculateScaledDeltas( cell->parameters[CellParams::AMR_DB] = dB; cell->parameters[CellParams::AMR_ALPHA1] = alpha; cell->parameters[CellParams::AMR_ALPHA2] = cell->parameters[CellParams::DX] * J / (Bperp + EPS); // Epsilon in denominator so we don't get infinities + cell->parameters[CellParams::P_ANISOTROPY] = Panisotropy; + // Experimental, current scaling is bulk velocity + cell->parameters[CellParams::AMR_VORTICITY] = amr_vorticity; } /*! \brief High-level scaled gradient calculation wrapper function. diff --git a/fieldsolver/gridGlue.cpp b/fieldsolver/gridGlue.cpp index a8073bdb0..4574d1e67 100644 --- a/fieldsolver/gridGlue.cpp +++ b/fieldsolver/gridGlue.cpp @@ -229,6 +229,7 @@ void getFieldsFromFsGrid( FsGrid< std::array, FS_STENCIL_WIDTH> & volumeFieldsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, dccrg::Dccrg& mpiGrid, const std::vector& cells @@ -311,7 +312,9 @@ void getFieldsFromFsGrid( std::array * volcell = volumeFieldsGrid.get(fsgridCell); std::array * bgcell = BgBGrid.get(fsgridCell); std::array * egradpecell = EGradPeGrid.get(fsgridCell); + std::array * dMomentscell = dMomentsGrid.get(fsgridCell); + // TODO consider pruning these and communicating only when required sendBuffer[ii].sums[FieldsToCommunicate::PERBXVOL] += volcell->at(fsgrids::volfields::PERBXVOL); sendBuffer[ii].sums[FieldsToCommunicate::PERBYVOL] += volcell->at(fsgrids::volfields::PERBYVOL); sendBuffer[ii].sums[FieldsToCommunicate::PERBZVOL] += volcell->at(fsgrids::volfields::PERBZVOL); @@ -324,6 +327,15 @@ void getFieldsFromFsGrid( sendBuffer[ii].sums[FieldsToCommunicate::dPERBZVOLdx] += volcell->at(fsgrids::volfields::dPERBZVOLdx) / technicalGrid.DX; sendBuffer[ii].sums[FieldsToCommunicate::dPERBZVOLdy] += volcell->at(fsgrids::volfields::dPERBZVOLdy) / technicalGrid.DY; sendBuffer[ii].sums[FieldsToCommunicate::dPERBZVOLdz] += volcell->at(fsgrids::volfields::dPERBZVOLdz) / technicalGrid.DZ; + sendBuffer[ii].sums[FieldsToCommunicate::dVxdx] += dMomentscell->at(fsgrids::dmoments::dVxdx) / technicalGrid.DX; + sendBuffer[ii].sums[FieldsToCommunicate::dVxdy] += dMomentscell->at(fsgrids::dmoments::dVxdy) / technicalGrid.DY; + sendBuffer[ii].sums[FieldsToCommunicate::dVxdz] += dMomentscell->at(fsgrids::dmoments::dVxdz) / technicalGrid.DZ; + sendBuffer[ii].sums[FieldsToCommunicate::dVydx] += dMomentscell->at(fsgrids::dmoments::dVydx) / technicalGrid.DX; + sendBuffer[ii].sums[FieldsToCommunicate::dVydy] += dMomentscell->at(fsgrids::dmoments::dVydy) / technicalGrid.DY; + sendBuffer[ii].sums[FieldsToCommunicate::dVydz] += dMomentscell->at(fsgrids::dmoments::dVydz) / technicalGrid.DZ; + sendBuffer[ii].sums[FieldsToCommunicate::dVzdx] += dMomentscell->at(fsgrids::dmoments::dVzdx) / technicalGrid.DX; + sendBuffer[ii].sums[FieldsToCommunicate::dVzdy] += dMomentscell->at(fsgrids::dmoments::dVzdy) / technicalGrid.DY; + sendBuffer[ii].sums[FieldsToCommunicate::dVzdz] += dMomentscell->at(fsgrids::dmoments::dVzdz) / technicalGrid.DZ; sendBuffer[ii].sums[FieldsToCommunicate::BGBXVOL] += bgcell->at(fsgrids::bgbfield::BGBXVOL); sendBuffer[ii].sums[FieldsToCommunicate::BGBYVOL] += bgcell->at(fsgrids::bgbfield::BGBYVOL); sendBuffer[ii].sums[FieldsToCommunicate::BGBZVOL] += bgcell->at(fsgrids::bgbfield::BGBZVOL); @@ -383,6 +395,15 @@ void getFieldsFromFsGrid( mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdx] = cellAggregate.second.sums[FieldsToCommunicate::dPERBZVOLdx] / cellAggregate.second.cells; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdy] = cellAggregate.second.sums[FieldsToCommunicate::dPERBZVOLdy] / cellAggregate.second.cells; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdz] = cellAggregate.second.sums[FieldsToCommunicate::dPERBZVOLdz] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdx] = cellAggregate.second.sums[FieldsToCommunicate::dVxdx] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdy] = cellAggregate.second.sums[FieldsToCommunicate::dVxdy] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdz] = cellAggregate.second.sums[FieldsToCommunicate::dVxdz] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydx] = cellAggregate.second.sums[FieldsToCommunicate::dVydx] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydy] = cellAggregate.second.sums[FieldsToCommunicate::dVydy] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydz] = cellAggregate.second.sums[FieldsToCommunicate::dVydz] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdx] = cellAggregate.second.sums[FieldsToCommunicate::dVzdx] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdy] = cellAggregate.second.sums[FieldsToCommunicate::dVzdy] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdz] = cellAggregate.second.sums[FieldsToCommunicate::dVzdz] / cellAggregate.second.cells; cellParams[CellParams::BGBXVOL] = cellAggregate.second.sums[FieldsToCommunicate::BGBXVOL] / cellAggregate.second.cells; cellParams[CellParams::BGBYVOL] = cellAggregate.second.sums[FieldsToCommunicate::BGBYVOL] / cellAggregate.second.cells; cellParams[CellParams::BGBZVOL] = cellAggregate.second.sums[FieldsToCommunicate::BGBZVOL] / cellAggregate.second.cells; @@ -410,6 +431,15 @@ void getFieldsFromFsGrid( mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdx] = 0; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdy] = 0; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdz] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdx] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdy] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdz] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydx] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydy] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydz] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdx] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdy] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdz] = 0; cellParams[CellParams::BGBXVOL] = 0; cellParams[CellParams::BGBYVOL] = 0; cellParams[CellParams::BGBZVOL] = 0; diff --git a/fieldsolver/gridGlue.hpp b/fieldsolver/gridGlue.hpp index 43a03d417..f58d40dbc 100644 --- a/fieldsolver/gridGlue.hpp +++ b/fieldsolver/gridGlue.hpp @@ -36,6 +36,15 @@ enum FieldsToCommunicate { CURVATUREX, CURVATUREY, CURVATUREZ, + dVxdx, + dVxdy, + dVxdz, + dVydx, + dVydy, + dVydz, + dVzdx, + dVzdy, + dVzdz, N_FIELDSTOCOMMUNICATE }; @@ -64,12 +73,14 @@ void feedMomentsIntoFsGrid(dccrg::Dccrg& * * This function assumes that proper grid coupling has been set up. */ -void getFieldsFromFsGrid(FsGrid< std::array, FS_STENCIL_WIDTH> & volumeFieldsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - dccrg::Dccrg& mpiGrid, - const std::vector& cells +void getFieldsFromFsGrid( + FsGrid< std::array, FS_STENCIL_WIDTH> & volumeFieldsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, + dccrg::Dccrg& mpiGrid, + const std::vector& cells ); /*! Copy background B fields and store them into DCCRG diff --git a/grid.cpp b/grid.cpp index 6f6f8f8c0..d0f284d90 100644 --- a/grid.cpp +++ b/grid.cpp @@ -95,6 +95,7 @@ void initializeGrids( FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, @@ -356,7 +357,7 @@ void initializeGrids( volGrid.updateGhostCells(); fsGridGhostTimer.stop(); phiprof::Timer getFieldsTimer {"getFieldsFromFsGrid"}; - getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, technicalGrid, mpiGrid, cells); + getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, dMomentsGrid, technicalGrid, mpiGrid, cells); getFieldsTimer.stop(); setBTimer.stop(); diff --git a/grid.h b/grid.h index 0b4c6628a..7f389a819 100644 --- a/grid.h +++ b/grid.h @@ -41,6 +41,7 @@ void initializeGrids( FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, diff --git a/parameters.cpp b/parameters.cpp index b6b5eb4bb..f360fde4a 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -167,6 +167,13 @@ Real P::alpha1CoarsenThreshold = -1.0; bool P::useAlpha2 = true; Real P::alpha2RefineThreshold = 0.5; Real P::alpha2CoarsenThreshold = -1.0; +bool P::useVorticity = false; +Real P::vorticityRefineThreshold = 0.5; +Real P::vorticityCoarsenThreshold = -1.0; +bool P::useAnisotropy = false; +Real P::anisotropyRefineThreshold = 0.5; +Real P::anisotropyCoarsenThreshold = -1.0; +int P::anisotropyMaxReflevel = 2; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; Real P::alphaDPSqWeight = 1.0; @@ -409,7 +416,7 @@ bool P::addParameters() { "ig_precipitation ig_deltaphi "+ "ig_inplanecurrent ig_b ig_e vg_drift vg_ionospherecoupling vg_connection vg_fluxrope fg_curvature "+ "vg_amr_drho vg_amr_du vg_amr_dpsq vg_amr_dbsq vg_amr_db vg_amr_alpha1 vg_amr_reflevel vg_amr_alpha2 "+ - "vg_amr_translate_comm vg_gridcoordinates fg_gridcoordinates "); + "vg_amr_translate_comm vg_gridcoordinates fg_gridcoordinates vg_pressure_anisotropy vg_amr_vorticity"); RP::addComposing( "variables_deprecated.output", @@ -484,6 +491,13 @@ bool P::addParameters() { RP::add("AMR.use_alpha2","Use J/B_perp as a refinement index", true); RP::add("AMR.alpha2_refine_threshold","Determines the minimum value of alpha_2 to refine cells", 0.5); RP::add("AMR.alpha2_coarsen_threshold","Determines the maximum value of alpha_2 to unrefine cells, default half of the refine threshold", -1.0); + RP::add("AMR.use_vorticity","Use vorticity as a refinement index", false); + RP::add("AMR.vorticity_refine_threshold","Determines the minimum value of vorticity to refine cells", 0.5); + RP::add("AMR.vorticity_coarsen_threshold","Determines the maximum value of vorticity to unrefine cells, default half of the refine threshold", -1.0); + RP::add("AMR.use_anisotropy","Use pressure anisotropy as a refinement index", false); + RP::add("AMR.anisotropy_refine_threshold","Determines the maximum value of pressure anisotropy to refine cells", 0.5); + RP::add("AMR.anisotropy_coarsen_threshold","Determines the minimum value of pressure anisotropy to unrefine cells, default twice the refine threshold", -1.0); + RP::add("AMR.anisotropy_max_reflevel","When anisotropy is below the refine threshold, defines the maximum level to refine to", 2); RP::add("AMR.refine_cadence","Refine every nth load balance", 5); RP::add("AMR.refine_after","Start refinement after this many simulation seconds", 0.0); RP::add("AMR.refine_radius","Maximum distance from origin to allow refinement within. Only induced refinement allowed outside this radius.", LARGE_REAL); @@ -788,6 +802,7 @@ void Parameters::getParameters() { } MPI_Abort(MPI_COMM_WORLD, 1); } + RP::get("AMR.use_alpha2",P::useAlpha2); RP::get("AMR.alpha2_refine_threshold",P::alpha2RefineThreshold); RP::get("AMR.alpha2_coarsen_threshold",P::alpha2CoarsenThreshold); @@ -801,6 +816,33 @@ void Parameters::getParameters() { MPI_Abort(MPI_COMM_WORLD, 1); } + RP::get("AMR.use_vorticity",P::useVorticity); + RP::get("AMR.vorticity_refine_threshold",P::vorticityRefineThreshold); + RP::get("AMR.vorticity_coarsen_threshold",P::vorticityCoarsenThreshold); + if (P::useVorticity && P::vorticityCoarsenThreshold < 0) { + P::vorticityCoarsenThreshold = P::vorticityRefineThreshold / 2.0; + } + if (P::useVorticity && P::vorticityRefineThreshold < 0) { + if (myRank == MASTER_RANK) { + cerr << "ERROR invalid vorticity refine threshold" << endl; + } + MPI_Abort(MPI_COMM_WORLD, 1); + } + + RP::get("AMR.use_anisotropy",P::useAnisotropy); + RP::get("AMR.anisotropy_refine_threshold", P::anisotropyRefineThreshold); + RP::get("AMR.anisotropy_coarsen_threshold", P::anisotropyCoarsenThreshold); + RP::get("AMR.anisotropy_max_reflevel", P::anisotropyMaxReflevel); + if (P::useAnisotropy && P::anisotropyCoarsenThreshold < 0) { + P::anisotropyCoarsenThreshold = P::anisotropyRefineThreshold * 2.0; + } + if (P::useAnisotropy && P::anisotropyRefineThreshold < 0) { + if (myRank == MASTER_RANK) { + cerr << "ERROR invalid anisotropy refine threshold" << endl; + } + MPI_Abort(MPI_COMM_WORLD, 1); + } + RP::get("AMR.refine_cadence",P::refineCadence); RP::get("AMR.refine_after",P::refineAfter); RP::get("AMR.refine_radius",P::refineRadius); diff --git a/parameters.h b/parameters.h index 87838ef9f..bd63bdbad 100644 --- a/parameters.h +++ b/parameters.h @@ -198,6 +198,14 @@ struct Parameters { static bool useAlpha2; static Real alpha2RefineThreshold; static Real alpha2CoarsenThreshold; + // TODO: consider renaming to alpha3 or something to that effect + static bool useVorticity; + static Real vorticityRefineThreshold; + static Real vorticityCoarsenThreshold; + static bool useAnisotropy; + static Real anisotropyRefineThreshold; + static Real anisotropyCoarsenThreshold; + static int anisotropyMaxReflevel; static uint refineCadence; static Real refineAfter; static Real refineRadius; diff --git a/projects/project.cpp b/projects/project.cpp index a4b174fec..801709ea6 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -563,13 +563,88 @@ namespace projects { return cell->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY && (cell->sysBoundaryLayer == 0 || cell->sysBoundaryLayer > 2); } + bool Project::shouldRefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const { + // Evaluate possible refinement for this cell + + // Cells too far from the ionosphere should not be refined but + // induced refinement still possible just beyond this r_max2 limit. + + std::array xyz {mpiGrid.get_center(id)}; + SpatialCell* cell {mpiGrid[id]}; + int refLevel {mpiGrid.get_refinement_level(id)}; + Real r2 {pow(xyz[0], 2) + pow(xyz[1], 2) + pow(xyz[2], 2)}; + + bool alpha1ShouldRefine = (P::useAlpha1 && cell->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold); + bool alpha2ShouldRefine = (P::useAlpha2 && cell->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold); + bool vorticityShouldRefine = (P::useVorticity && cell->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold); + bool anisotropyShouldRefine = (P::useAnisotropy && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel); + + bool shouldRefine { + (r2 < r_max2) && ( + alpha1ShouldRefine || + alpha2ShouldRefine || + vorticityShouldRefine || + anisotropyShouldRefine + ) + }; + + if( + // If this cell is planned to be refined, but is outside the allowed refinement region, cancel that refinement. + // Induced refinement still possible just beyond that limit. + (xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) + || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) + || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ)) { + shouldRefine = false; + } + + return shouldRefine; + } + + bool Project::shouldUnrefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const { + // Evaluate possible unrefinement for this cell + + // Cells too far from the ionosphere should be unrefined but + // induced refinement still possible just beyond this r_max2 limit. + + std::array xyz {mpiGrid.get_center(id)}; + SpatialCell* cell {mpiGrid[id]}; + int refLevel {mpiGrid.get_refinement_level(id)}; + Real r2 {pow(xyz[0], 2) + pow(xyz[1], 2) + pow(xyz[2], 2)}; + + bool alpha1ShouldUnrefine = (!P::useAlpha1 || cell->parameters[CellParams::AMR_ALPHA1] < P::alpha1CoarsenThreshold); + bool alpha2ShouldUnrefine = (!P::useAlpha2 || cell->parameters[CellParams::AMR_ALPHA2] < P::alpha2CoarsenThreshold); + bool vorticityShouldUnrefine = (!P::useVorticity || cell->parameters[CellParams::AMR_VORTICITY] < P::vorticityCoarsenThreshold); + bool anisotropyShouldUnrefine = (!P::useAnisotropy || cell->parameters[CellParams::P_ANISOTROPY] > P::anisotropyCoarsenThreshold || refLevel > P::anisotropyMaxReflevel); + + bool shouldUnrefine { + (r2 > r_max2) || ( + alpha1ShouldUnrefine && + alpha2ShouldUnrefine && + vorticityShouldUnrefine && + anisotropyShouldUnrefine + ) + }; + + if( + // If this cell is planned to remain at the current refinement level, but is outside the allowed refinement region, + // attempt to unrefine it instead. (If it is already at the lowest refinement level, DCCRG should not go belly-up.) + // Induced refinement still possible just beyond that limit. + (xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) + || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) + || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ)) { + shouldUnrefine = true; + } + + return shouldUnrefine; + } + int Project::adaptRefinement( dccrg::Dccrg& mpiGrid ) const { phiprof::Timer refinesTimer {"Set refines"}; int myRank; MPI_Comm_rank(MPI_COMM_WORLD,&myRank); int refines {0}; - if (!P::useAlpha1 && !P::useAlpha2) { + if (!P::useAlpha1 && !P::useAlpha2 && !P::useAnisotropy && !P::useVorticity) { if (myRank == MASTER_RANK) { std::cout << "WARNING All refinement indices disabled" << std::endl; } @@ -581,10 +656,7 @@ namespace projects { //#pragma omp parallel for for (CellID id : cells) { - std::array xyz {mpiGrid.get_center(id)}; - SpatialCell* cell {mpiGrid[id]}; int refLevel {mpiGrid.get_refinement_level(id)}; - Real r2 {pow(xyz[0], 2) + pow(xyz[1], 2) + pow(xyz[2], 2)}; if (!canRefine(mpiGrid[id])) { // Skip refining, touching boundaries during runtime breaks everything @@ -593,58 +665,18 @@ namespace projects { } else { // Evaluate possible refinement or unrefinement for this cell - // Cells too far from the ionosphere should be unrefined but - // induced refinement still possible just beyond this r_max2 limit. - bool shouldRefine {(r2 < r_max2) && ((P::useAlpha1 ? cell->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold : false) || (P::useAlpha2 ? cell->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold : false))}; - bool shouldUnrefine {(r2 > r_max2) || ((P::useAlpha1 ? cell->parameters[CellParams::AMR_ALPHA1] < P::alpha1CoarsenThreshold : true) && (P::useAlpha2 ? cell->parameters[CellParams::AMR_ALPHA2] < P::alpha2CoarsenThreshold : true))}; - - if(shouldRefine - // If this cell is planned to be refined, but is outside the allowed refinement region, cancel that refinement. - // Induced refinement still possible just beyond that limit. - && ((xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) - || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) - || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ))) { - shouldRefine = false; - } - if(!shouldUnrefine - // If this cell is planned to remain at the current refinement level, but is outside the allowed refinement region, - // attempt to unrefine it instead. (If it is already at the lowest refinement level, DCCRG should not go belly-up.) - // Induced refinement still possible just beyond that limit. - && ((xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) - || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) - || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ))) { - shouldUnrefine = true; - } + bool shouldRefine = shouldRefineCell(mpiGrid, id, r_max2); + bool shouldUnrefine = shouldUnrefineCell(mpiGrid, id, r_max2); // Finally, check neighbors int refined_neighbors {0}; int coarser_neighbors {0}; for (const auto& [neighbor, dir] : mpiGrid.get_face_neighbors_of(id)) { // Evaluate all face neighbors of the current cell - std::array neighborXyz {mpiGrid.get_center(neighbor)}; - Real neighborR2 {pow(neighborXyz[0], 2) + pow(neighborXyz[1], 2) + pow(neighborXyz[2], 2)}; - const int neighborRef = mpiGrid.get_refinement_level(neighbor); - // Induced refinement still possible just beyond the r_max2 limit. - bool shouldRefineNeighbor {(neighborR2 < r_max2) && ((P::useAlpha1 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold : false) || (P::useAlpha2 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold : false))}; - if(shouldRefineNeighbor && - // If the neighbor is planned to be refined, but is outside the allowed refinement region, cancel that refinement. - // Induced refinement still possible just beyond that limit. - ((neighborXyz[0] < P::refinementMinX) || (neighborXyz[0] > P::refinementMaxX) - || (neighborXyz[1] < P::refinementMinY) || (neighborXyz[1] > P::refinementMaxY) - || (neighborXyz[2] < P::refinementMinZ) || (neighborXyz[2] > P::refinementMaxZ))) { - shouldRefineNeighbor = false; - } - // Induced refinement still possible just beyond the r_max2 limit. - bool shouldUnrefineNeighbor {(neighborR2 > r_max2) || ((P::useAlpha1 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA1] < P::alpha1CoarsenThreshold : true) && (P::useAlpha2 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA2] < P::alpha2CoarsenThreshold : true))}; - if(!shouldUnrefineNeighbor && - // If the neighbor is planned to remain at the current refinement level, but is outside the allowed refinement region, - // consider it as unrefining instead for purposes of evaluating the neighbors of this cell. - // Induced refinement still possible just beyond that limit. - ((neighborXyz[0] < P::refinementMinX) || (neighborXyz[0] > P::refinementMaxX) - || (neighborXyz[1] < P::refinementMinY) || (neighborXyz[1] > P::refinementMaxY) - || (neighborXyz[2] < P::refinementMinZ) || (neighborXyz[2] > P::refinementMaxZ))) { - shouldUnrefineNeighbor = true; - } + bool shouldRefineNeighbor = shouldRefineCell(mpiGrid, neighbor, r_max2); + bool shouldUnrefineNeighbor = shouldUnrefineCell(mpiGrid, neighbor, r_max2); + int neighborRef {mpiGrid.get_refinement_level(neighbor)}; + if (neighborRef > refLevel && !shouldUnrefineNeighbor) { ++refined_neighbors; } else if (neighborRef < refLevel && !shouldRefineNeighbor) { diff --git a/projects/project.h b/projects/project.h index 0e7f44dc3..ecf1fa38f 100644 --- a/projects/project.h +++ b/projects/project.h @@ -85,6 +85,10 @@ namespace projects { virtual bool canRefine(spatial_cell::SpatialCell* cell) const; + virtual bool shouldRefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const; + + virtual bool shouldUnrefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const; + virtual bool refineSpatialCells( dccrg::Dccrg& mpiGrid ) const; /*!\brief Adapts refinement by one level according to the project. Returns true if any cells were refined, false if not. diff --git a/spatial_cell_cpu.cpp b/spatial_cell_cpu.cpp index b7f0f2c0d..3fcb17d5d 100644 --- a/spatial_cell_cpu.cpp +++ b/spatial_cell_cpu.cpp @@ -743,7 +743,7 @@ namespace spatial_cell { // Refinement parameters if ((SpatialCell::mpi_transfer_type & Transfer::REFINEMENT_PARAMETERS)){ displacements.push_back(reinterpret_cast(this->parameters.data() + CellParams::AMR_ALPHA1) - reinterpret_cast(this)); - block_lengths.push_back(sizeof(Real) * (CellParams::AMR_ALPHA2 - CellParams::AMR_ALPHA1 + 1)); // This is just 2, but let's be explicit + block_lengths.push_back(sizeof(Real) * (CellParams::AMR_VORTICITY - CellParams::AMR_ALPHA1 + 1)); // This is just 4, but let's be explicit. } // Copy random number generator state variables diff --git a/spatial_cell_cpu.hpp b/spatial_cell_cpu.hpp index 55281d2fc..d9ec4f5ee 100644 --- a/spatial_cell_cpu.hpp +++ b/spatial_cell_cpu.hpp @@ -316,6 +316,7 @@ namespace spatial_cell { //random_data* get_rng_data_buffer(); // Member variables // + std::array derivativesV; // Derivatives of V for AMR std::array derivativesBVOL; /**< Derivatives of BVOL needed by the acceleration. * Separate array because it does not need to be communicated.*/ //Real parameters[CellParams::N_SPATIAL_CELL_PARAMS]; /**< Bulk variables in this spatial cell.*/ diff --git a/vlasiator.cpp b/vlasiator.cpp index c0af9d8a3..1e6edc9a6 100644 --- a/vlasiator.cpp +++ b/vlasiator.cpp @@ -481,6 +481,7 @@ int simulate(int argn,char* args[]) { BgBGrid, momentsGrid, momentsDt2Grid, + dMomentsGrid, EGrid, EGradPeGrid, volGrid, @@ -612,7 +613,7 @@ int simulate(int argn,char* args[]) { phiprof::Timer getFieldsTimer {"getFieldsFromFsGrid"}; volGrid.updateGhostCells(); - getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, technicalGrid, mpiGrid, cells); + getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, dMomentsGrid, technicalGrid, mpiGrid, cells); getFieldsTimer.stop(); // Build communicator for ionosphere solving @@ -751,7 +752,10 @@ int simulate(int argn,char* args[]) { CellParams::RHOQ, CellParams::P_11, CellParams::P_22, - CellParams::P_33 + CellParams::P_33, + CellParams::P_23, + CellParams::P_13, + CellParams::P_12 ); computeMomentsTimer.stop(); } @@ -1166,7 +1170,10 @@ int simulate(int argn,char* args[]) { CellParams::RHOQ_DT2, CellParams::P_11_DT2, CellParams::P_22_DT2, - CellParams::P_33_DT2 + CellParams::P_33_DT2, + CellParams::P_23_DT2, + CellParams::P_13_DT2, + CellParams::P_12_DT2 ); momentsTimer.stop(); @@ -1207,7 +1214,7 @@ int simulate(int argn,char* args[]) { // Copy results back from fsgrid. volGrid.updateGhostCells(); technicalGrid.updateGhostCells(); - getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, technicalGrid, mpiGrid, cells); + getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, dMomentsGrid, technicalGrid, mpiGrid, cells); getFieldsTimer.stop(); propagateTimer.stop(cells.size(),"SpatialCells"); addTimedBarrier("barrier-after-field-solver"); @@ -1276,7 +1283,10 @@ int simulate(int argn,char* args[]) { CellParams::RHOQ, CellParams::P_11, CellParams::P_22, - CellParams::P_33 + CellParams::P_33, + CellParams::P_23, + CellParams::P_13, + CellParams::P_12 ); momentsTimer.stop(); diff --git a/vlasovmover.h b/vlasovmover.h index 9289d8396..85dd72401 100644 --- a/vlasovmover.h +++ b/vlasovmover.h @@ -64,7 +64,10 @@ void calculateInterpolatedVelocityMoments( const int cp_rhoq, const int cp_p11, const int cp_p22, - const int cp_p33 + const int cp_p33, + const int cp_p23, + const int cp_p13, + const int cp_p12 ); /*! diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 995dc12f5..7e032674a 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -56,6 +56,9 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, cell->parameters[CellParams::P_11] = 0.0; cell->parameters[CellParams::P_22] = 0.0; cell->parameters[CellParams::P_33] = 0.0; + cell->parameters[CellParams::P_23] = 0.0; + cell->parameters[CellParams::P_13] = 0.0; + cell->parameters[CellParams::P_12] = 0.0; } // Loop over all particle species @@ -125,10 +128,7 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, const Real mass = getObjectWrapper().particleSpecies[popID].mass; // Temporary array for storing moments - Real array[3]; - for (int i=0; i<3; ++i) { - array[i] = 0.0; - } + std::vector array(6, 0.0); // Calculate species' contribution to second velocity moments Population & pop = cell->get_population(popID); @@ -147,9 +147,12 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, pop.P[2] = mass*array[2]; if (!computePopulationMomentsOnly) { - cell->parameters[CellParams::P_11] += pop.P[0]; - cell->parameters[CellParams::P_22] += pop.P[1]; - cell->parameters[CellParams::P_33] += pop.P[2]; + cell->parameters[CellParams::P_11] += mass * array[0]; + cell->parameters[CellParams::P_22] += mass * array[1]; + cell->parameters[CellParams::P_33] += mass * array[2]; + cell->parameters[CellParams::P_23] += mass * array[3]; + cell->parameters[CellParams::P_13] += mass * array[4]; + cell->parameters[CellParams::P_12] += mass * array[5]; } } // for-loop over particle species } @@ -187,6 +190,9 @@ void calculateMoments_R( cell->parameters[CellParams::P_11_R] = 0.0; cell->parameters[CellParams::P_22_R] = 0.0; cell->parameters[CellParams::P_33_R] = 0.0; + cell->parameters[CellParams::P_23_R] = 0.0; + cell->parameters[CellParams::P_13_R] = 0.0; + cell->parameters[CellParams::P_12_R] = 0.0; } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); @@ -277,10 +283,7 @@ void calculateMoments_R( const Real mass = getObjectWrapper().particleSpecies[popID].mass; // Temporary array where species' contribution to 2nd moments is accumulated - Real array[3]; - for (int i=0; i<3; ++i) { - array[i] = 0.0; - } + std::vector array(6, 0.0); // Calculate species' contribution to second velocity moments Population & pop = cell->get_population(popID); @@ -298,9 +301,12 @@ void calculateMoments_R( pop.P_R[1] = mass*array[1]; pop.P_R[2] = mass*array[2]; - cell->parameters[CellParams::P_11_R] += pop.P_R[0]; - cell->parameters[CellParams::P_22_R] += pop.P_R[1]; - cell->parameters[CellParams::P_33_R] += pop.P_R[2]; + cell->parameters[CellParams::P_11_R] += mass * array[0]; + cell->parameters[CellParams::P_22_R] += mass * array[1]; + cell->parameters[CellParams::P_33_R] += mass * array[2]; + cell->parameters[CellParams::P_23_R] += mass * array[3]; + cell->parameters[CellParams::P_13_R] += mass * array[4]; + cell->parameters[CellParams::P_12_R] += mass * array[5]; } // for-loop over spatial cells } // for-loop over particle species @@ -341,6 +347,9 @@ void calculateMoments_V( cell->parameters[CellParams::P_11_V] = 0.0; cell->parameters[CellParams::P_22_V] = 0.0; cell->parameters[CellParams::P_33_V] = 0.0; + cell->parameters[CellParams::P_23_V] = 0.0; + cell->parameters[CellParams::P_13_V] = 0.0; + cell->parameters[CellParams::P_12_V] = 0.0; } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); @@ -417,10 +426,7 @@ void calculateMoments_V( const Real mass = getObjectWrapper().particleSpecies[popID].mass; // Temporary array where moments are stored - Real array[3]; - for (int i=0; i<3; ++i) { - array[i] = 0.0; - } + std::vector array(6, 0.0); // Calculate species' contribution to second velocity moments Population & pop = cell->get_population(popID); @@ -439,9 +445,12 @@ void calculateMoments_V( pop.P_V[1] = mass*array[1]; pop.P_V[2] = mass*array[2]; - cell->parameters[CellParams::P_11_V] += pop.P_V[0]; - cell->parameters[CellParams::P_22_V] += pop.P_V[1]; - cell->parameters[CellParams::P_33_V] += pop.P_V[2]; + cell->parameters[CellParams::P_11_V] += mass * array[0]; + cell->parameters[CellParams::P_22_V] += mass * array[1]; + cell->parameters[CellParams::P_33_V] += mass * array[2]; + cell->parameters[CellParams::P_23_V] += mass * array[3]; + cell->parameters[CellParams::P_13_V] += mass * array[4]; + cell->parameters[CellParams::P_12_V] += mass * array[5]; } // for-loop over spatial cells } // for-loop over particle species diff --git a/vlasovsolver/cpu_moments.h b/vlasovsolver/cpu_moments.h index ce98931cb..962920c41 100644 --- a/vlasovsolver/cpu_moments.h +++ b/vlasovsolver/cpu_moments.h @@ -43,7 +43,7 @@ void blockVelocityFirstMoments(const Realf* avgs,const Real* blockParams, template void blockVelocitySecondMoments(const Realf* avgs,const Real* blockParams, const REAL v[3], - REAL* array); + std::vector& array); void calculateMoments_R(dccrg::Dccrg& mpiGrid, const std::vector& cells, @@ -114,13 +114,16 @@ void blockVelocitySecondMoments( const REAL averageVX, const REAL averageVY, const REAL averageVZ, - REAL* array) { + std::vector& array) { const Real HALF = 0.5; Real nvx2_sum = 0.0; Real nvy2_sum = 0.0; Real nvz2_sum = 0.0; + Real nvyvz_sum = 0.0; + Real nvxvz_sum = 0.0; + Real nvxvy_sum = 0.0; for (uint k=0; k& cells = getLocalCells(); @@ -537,6 +540,9 @@ void calculateInterpolatedVelocityMoments( SC->parameters[cp_p11] = 0.5* ( SC->parameters[CellParams::P_11_R] + SC->parameters[CellParams::P_11_V] ); SC->parameters[cp_p22] = 0.5* ( SC->parameters[CellParams::P_22_R] + SC->parameters[CellParams::P_22_V] ); SC->parameters[cp_p33] = 0.5* ( SC->parameters[CellParams::P_33_R] + SC->parameters[CellParams::P_33_V] ); + SC->parameters[cp_p23] = 0.5* ( SC->parameters[CellParams::P_23_R] + SC->parameters[CellParams::P_23_V] ); + SC->parameters[cp_p13] = 0.5* ( SC->parameters[CellParams::P_13_R] + SC->parameters[CellParams::P_13_V] ); + SC->parameters[cp_p12] = 0.5* ( SC->parameters[CellParams::P_12_R] + SC->parameters[CellParams::P_12_V] ); for (uint popID=0; popIDget_population(popID);