Skip to content

Relative entropy #1329

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

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

Relative entropy #1329

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
775745a
Add Boltzman H function
Jul 27, 2026
bb463f2
Add relative entropy
Jul 27, 2026
825f2bf
Relative entropy per particle
Jul 27, 2026
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
36 changes: 36 additions & 0 deletions datareduction/datareducer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -825,6 +825,42 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti
continue;
}
}
if(P::systemWriteAllDROs || lowercase == "populations_entropy" || lowercase == "populations_vg_entropy") {
// Per-population Boltzmann H function ~ f log f d3v
for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) {
species::Species& species=getObjectWrapper().particleSpecies[i];
const std::string& pop = species.name;
outputReducer->addOperator(new DRO::VariableEntropy(i));
outputReducer->addMetadata(outputReducer->size()-1,"m^-3","\\mathrm{m}^{-1}","$H_\\mathrm{"+pop+"}$","1.0");
}
if(!P::systemWriteAllDROs) {
continue;
}
}
if(P::systemWriteAllDROs || lowercase == "populations_relative_entropy" || lowercase == "populations_vg_relative_entropy") {
// Per-population relative entropy ~ f log (f/f_M) d3v, where f_M is a Maxwellian with same density n, flow speed U and three component temperature T
for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) {
species::Species& species=getObjectWrapper().particleSpecies[i];
const std::string& pop = species.name;
outputReducer->addOperator(new DRO::VariableRelativeEntropy(i));
outputReducer->addMetadata(outputReducer->size()-1,"m^-3","\\mathrm{m}^{-1}","$\\sigma_\\mathrm{"+pop+",rel}$","1.0");
}
if(!P::systemWriteAllDROs) {
continue;
}
}
if(P::systemWriteAllDROs || lowercase == "populations_nonmaxwellianess" || lowercase == "populations_vg_nonmaxwellianess") {
// Per-population relative entropy per particle ~ f/n log (f/f_M) d3v, where f_M is a Maxwellian with same density n, flow speed U and three component temperature T
for(unsigned int i =0; i < getObjectWrapper().particleSpecies.size(); i++) {
species::Species& species=getObjectWrapper().particleSpecies[i];
const std::string& pop = species.name;
outputReducer->addOperator(new DRO::VariableNonMaxwellianess(i));
outputReducer->addMetadata(outputReducer->size()-1,"","","$\\sigma_\\mathrm{"+pop+",rel} \\,/\\, n_{"+pop+"}$","1.0");
}
if(!P::systemWriteAllDROs) {
continue;
}
}
if(P::systemWriteAllDROs || lowercase == "bvolderivs" || lowercase == "b_vol_derivs" || lowercase == "b_vol_derivatives" || lowercase == "vg_b_vol_derivatives" || lowercase == "derivs") {
// Volume-averaged derivatives
outputReducer->addOperator(new DRO::DataReductionOperatorBVOLDerivatives("vg_derivatives/vg_dperbxvoldx",bvolderivatives::dPERBXVOLdx,1));
Expand Down
291 changes: 284 additions & 7 deletions datareduction/datareductionoperator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -758,10 +758,287 @@ namespace DRO {
return true;
}

/********
Next level of helper functions - these include threading and calculate zeroth or first velocity moments or the
diagonal / off-diagonal pressure tensor components for
thermal or non-thermal populations ********/
VariableEntropy::VariableEntropy(cuint _popID): DataReductionOperator(),popID(_popID) {
popName = getObjectWrapper().particleSpecies[popID].name;
}
VariableEntropy::~VariableEntropy() { }

std::string VariableEntropy::getName() const {return popName + "/vg_entropy";}

bool VariableEntropy::getDataVectorInfo(std::string& dataType,unsigned int& dataSize,unsigned int& vectorSize) const {
dataType = "float";
dataSize = sizeof(Real);
vectorSize = 1;
return true;
}

bool VariableEntropy::reduceData(const SpatialCell* cell,char* buffer) {
#ifdef USE_GPU
const vmesh::VelocityBlockContainer* VBC = cell->dev_get_velocity_blocks(popID);
#else
const vmesh::VelocityBlockContainer* VBC = cell->get_velocity_blocks(popID);
#endif
// ARCH interface includes OpenMP looping including critical regions for thread summation
{
Real sum = 0.0;

if (cell->get_number_of_velocity_blocks(popID) != 0)
arch::parallel_reduce<arch::sum>({WID, WID, WID, (uint)cell->get_number_of_velocity_blocks(popID)},
ARCH_LOOP_LAMBDA(const uint i, const uint j, const uint k, const uint n, Real *lsum ) {

const Realf *block_data = VBC->getData(n);
const Real *block_parameters = VBC->getParameters(n);
const Real DV3 = block_parameters[BlockParams::DVX]
* block_parameters[BlockParams::DVY] * block_parameters[BlockParams::DVZ];

Real inc = 0.0;
if (block_data[cellIndex(i,j,k)] > 0.) {
inc = block_data[cellIndex(i,j,k)] * log(block_data[cellIndex(i,j,k)]) * DV3;
}
lsum[0] += inc;
}, sum);

Entropy = sum;
}
const char* ptr = reinterpret_cast<const char*>(&Entropy);
for (uint i = 0; i < sizeof(Real); ++i) buffer[i] = ptr[i];
return true;
}

bool VariableEntropy::setSpatialCell(const SpatialCell* cell) {
Entropy = 0.0;
return true;
}

/*! \brief Relative Entropy
* Calculate for a population the relative entropy as defined by Cassak et al. *2023) as
* sigma_M = integral f log (f/f_M) d3v
* where
* f: measured VDF for a population
* f_M: a Maxwellian with same density n, flow speed U and three component temperature T
* No extra parameters need to be set in cfg file.
*/
VariableRelativeEntropy::VariableRelativeEntropy(cuint _popID): DataReductionOperator(),popID(_popID) {
popName = getObjectWrapper().particleSpecies[popID].name;
speciesMass = getObjectWrapper().particleSpecies[popID].mass;
}
VariableRelativeEntropy::~VariableRelativeEntropy() { }

std::string VariableRelativeEntropy::getName() const {return popName + "/vg_relative_entropy";}

bool VariableRelativeEntropy::getDataVectorInfo(std::string& dataType,unsigned int& dataSize,unsigned int& vectorSize) const {
dataType = "float";
dataSize = sizeof(Real);
vectorSize = 1;
return true;
}

bool VariableRelativeEntropy::reduceData(const SpatialCell* cell,char* buffer) {
#ifdef USE_GPU
const vmesh::VelocityBlockContainer* VBC = cell->dev_get_velocity_blocks(popID);
#else
const vmesh::VelocityBlockContainer* VBC = cell->get_velocity_blocks(popID);
#endif

// First pass we accumulate diagonal terms of the pressure tensor and divide by density to get temperatures
const Real HALF = 0.5;
{
Real sum[3] = {0.0, 0.0, 0.0};
Real averageVX = this->averageVX, averageVY = this->averageVY, averageVZ = this->averageVZ;

if (cell->get_number_of_velocity_blocks(popID) != 0)
arch::parallel_reduce<arch::sum>({WID, WID, WID, (uint)cell->get_number_of_velocity_blocks(popID)},
ARCH_LOOP_LAMBDA(const uint i, const uint j, const uint k, const uint n, Real *lsum ){

const Realf *block_data = VBC->getData(n);
const Real *block_parameters = VBC->getParameters(n);
const Real VX = block_parameters[BlockParams::VXCRD] + (i + HALF)*block_parameters[BlockParams::DVX];
const Real VY = block_parameters[BlockParams::VYCRD] + (j + HALF)*block_parameters[BlockParams::DVY];
const Real VZ = block_parameters[BlockParams::VZCRD] + (k + HALF)*block_parameters[BlockParams::DVZ];
const Real DV3 = block_parameters[BlockParams::DVX]
* block_parameters[BlockParams::DVY] * block_parameters[BlockParams::DVZ];

lsum[0] += block_data[cellIndex(i,j,k)] * (VX - averageVX) * (VX - averageVX) * DV3;
lsum[1] += block_data[cellIndex(i,j,k)] * (VY - averageVY) * (VY - averageVY) * DV3;
lsum[2] += block_data[cellIndex(i,j,k)] * (VZ - averageVZ) * (VZ - averageVZ) * DV3;
}, sum);

Temperature[0] = sum[0] / density; // this really is k_B T / m
Temperature[1] = sum[1] / density;
Temperature[2] = sum[2] / density;
}
// Second pass we compute the relative entropy
{
Real sum = 0.0;
Real density = this->density;
Real averageVX = this->averageVX, averageVY = this->averageVY, averageVZ = this->averageVZ;
Real* Temperature = this->Temperature;

if (cell->get_number_of_velocity_blocks(popID) != 0)
arch::parallel_reduce<arch::sum>({WID, WID, WID, (uint)cell->get_number_of_velocity_blocks(popID)},
ARCH_LOOP_LAMBDA(const uint i, const uint j, const uint k, const uint n, Real *lsum ) {

const Realf *block_data = VBC->getData(n);
const Real *block_parameters = VBC->getParameters(n);
const Real VX = block_parameters[BlockParams::VXCRD] + (i + HALF)*block_parameters[BlockParams::DVX];
const Real VY = block_parameters[BlockParams::VYCRD] + (j + HALF)*block_parameters[BlockParams::DVY];
const Real VZ = block_parameters[BlockParams::VZCRD] + (k + HALF)*block_parameters[BlockParams::DVZ];
const Real DV3 = block_parameters[BlockParams::DVX]
* block_parameters[BlockParams::DVY] * block_parameters[BlockParams::DVZ];

Real theta_x = (VX-averageVX) * (VX-averageVX) / (2.*Temperature[0]);
Real theta_y = (VY-averageVY) * (VY-averageVY) / (2.*Temperature[1]);
Real theta_z = (VZ-averageVZ) * (VZ-averageVZ) / (2.*Temperature[2]);
Real fM = density *
sqrt(1./(2.*M_PI*Temperature[0])) *
sqrt(1./(2.*M_PI*Temperature[1])) *
sqrt(1./(2.*M_PI*Temperature[2])) *
exp(-theta_x - theta_y - theta_z);

Real inc = 0.0;
if (block_data[cellIndex(i,j,k)] > 0.) {
inc = block_data[cellIndex(i,j,k)] * log(block_data[cellIndex(i,j,k)]/fM) * DV3;
}
lsum[0] += inc;
}, sum);

RelativeEntropy = sum;
}

const char* ptr = reinterpret_cast<const char*>(&RelativeEntropy);
for (uint i = 0; i < sizeof(Real); ++i) buffer[i] = ptr[i];
return true;
}

bool VariableRelativeEntropy::setSpatialCell(const SpatialCell* cell) {
density = cell-> parameters[CellParams::RHOM] / speciesMass;
averageVX = cell-> parameters[CellParams::VX];
averageVY = cell-> parameters[CellParams::VY];
averageVZ = cell-> parameters[CellParams::VZ];
for(int i = 0; i < 3; i++) {
Temperature[i] = 0.0;
}
RelativeEntropy = 0.0;
return true;
}

/*! \brief Non-Maxwellianess
* Calculate for a population the relative entropy per density as defined by Cassak et al. *2023) as
* sigma_M = integral f log (f/f_M) d3v / n
* where
* f: measured VDF for a population
* f_M: a Maxwellian with same density n, flow speed U and three component temperature T
* No extra parameters need to be set in cfg file.
*/
VariableNonMaxwellianess::VariableNonMaxwellianess(cuint _popID): DataReductionOperator(),popID(_popID) {
popName = getObjectWrapper().particleSpecies[popID].name;
speciesMass = getObjectWrapper().particleSpecies[popID].mass;
}
VariableNonMaxwellianess::~VariableNonMaxwellianess() { }

std::string VariableNonMaxwellianess::getName() const {return popName + "/vg_nonmaxwellianess";}

bool VariableNonMaxwellianess::getDataVectorInfo(std::string& dataType,unsigned int& dataSize,unsigned int& vectorSize) const {
dataType = "float";
dataSize = sizeof(Real);
vectorSize = 1;
return true;
}

bool VariableNonMaxwellianess::reduceData(const SpatialCell* cell,char* buffer) {
#ifdef USE_GPU
const vmesh::VelocityBlockContainer* VBC = cell->dev_get_velocity_blocks(popID);
#else
const vmesh::VelocityBlockContainer* VBC = cell->get_velocity_blocks(popID);
#endif

// First pass we accumulate diagonal terms of the pressure tensor and divide by density to get temperatures
const Real HALF = 0.5;
{
Real sum[3] = {0.0, 0.0, 0.0};
Real averageVX = this->averageVX, averageVY = this->averageVY, averageVZ = this->averageVZ;

if (cell->get_number_of_velocity_blocks(popID) != 0)
arch::parallel_reduce<arch::sum>({WID, WID, WID, (uint)cell->get_number_of_velocity_blocks(popID)},
ARCH_LOOP_LAMBDA(const uint i, const uint j, const uint k, const uint n, Real *lsum ){

const Realf *block_data = VBC->getData(n);
const Real *block_parameters = VBC->getParameters(n);
const Real VX = block_parameters[BlockParams::VXCRD] + (i + HALF)*block_parameters[BlockParams::DVX];
const Real VY = block_parameters[BlockParams::VYCRD] + (j + HALF)*block_parameters[BlockParams::DVY];
const Real VZ = block_parameters[BlockParams::VZCRD] + (k + HALF)*block_parameters[BlockParams::DVZ];
const Real DV3 = block_parameters[BlockParams::DVX]
* block_parameters[BlockParams::DVY] * block_parameters[BlockParams::DVZ];

lsum[0] += block_data[cellIndex(i,j,k)] * (VX - averageVX) * (VX - averageVX) * DV3;
lsum[1] += block_data[cellIndex(i,j,k)] * (VY - averageVY) * (VY - averageVY) * DV3;
lsum[2] += block_data[cellIndex(i,j,k)] * (VZ - averageVZ) * (VZ - averageVZ) * DV3;
}, sum);

Temperature[0] = sum[0] / density; // this really is k_B T / m
Temperature[1] = sum[1] / density;
Temperature[2] = sum[2] / density;
}
// Second pass we compute the relative entropy
{
Real sum = 0.0;
Real density = this->density;
Real averageVX = this->averageVX, averageVY = this->averageVY, averageVZ = this->averageVZ;
Real* Temperature = this->Temperature;

if (cell->get_number_of_velocity_blocks(popID) != 0)
arch::parallel_reduce<arch::sum>({WID, WID, WID, (uint)cell->get_number_of_velocity_blocks(popID)},
ARCH_LOOP_LAMBDA(const uint i, const uint j, const uint k, const uint n, Real *lsum ) {

const Realf *block_data = VBC->getData(n);
const Real *block_parameters = VBC->getParameters(n);
const Real VX = block_parameters[BlockParams::VXCRD] + (i + HALF)*block_parameters[BlockParams::DVX];
const Real VY = block_parameters[BlockParams::VYCRD] + (j + HALF)*block_parameters[BlockParams::DVY];
const Real VZ = block_parameters[BlockParams::VZCRD] + (k + HALF)*block_parameters[BlockParams::DVZ];
const Real DV3 = block_parameters[BlockParams::DVX]
* block_parameters[BlockParams::DVY] * block_parameters[BlockParams::DVZ];

Real theta_x = (VX-averageVX) * (VX-averageVX) / (2.*Temperature[0]);
Real theta_y = (VY-averageVY) * (VY-averageVY) / (2.*Temperature[1]);
Real theta_z = (VZ-averageVZ) * (VZ-averageVZ) / (2.*Temperature[2]);
Real fM = density *
sqrt(1./(2.*M_PI*Temperature[0])) *
sqrt(1./(2.*M_PI*Temperature[1])) *
sqrt(1./(2.*M_PI*Temperature[2])) *
exp(-theta_x - theta_y - theta_z);

Real inc = 0.0;
if (block_data[cellIndex(i,j,k)] > 0.) {
inc = block_data[cellIndex(i,j,k)] * log(block_data[cellIndex(i,j,k)]/fM) * DV3;
}
lsum[0] += inc;
}, sum);

// Divide out the density to get a dimensionless per particle quantity
NonMaxwellianess = divideIfNonZero(sum, density);
}

const char* ptr = reinterpret_cast<const char*>(&NonMaxwellianess);
for (uint i = 0; i < sizeof(Real); ++i) buffer[i] = ptr[i];
return true;
}

bool VariableNonMaxwellianess::setSpatialCell(const SpatialCell* cell) {
density = cell-> parameters[CellParams::RHOM] / speciesMass;
averageVX = cell-> parameters[CellParams::VX];
averageVY = cell-> parameters[CellParams::VY];
averageVZ = cell-> parameters[CellParams::VZ];
for(int i = 0; i < 3; i++) {
Temperature[i] = 0.0;
}
NonMaxwellianess = 0.0;
return true;
}

/********
Next level of helper functions - these include threading and calculate zeroth or first velocity moments or the
diagonal / off-diagonal pressure tensor components for
thermal or non-thermal populations ********/

//Calculates rho thermal or rho non-thermal
static void rhoNonthermalCalculation( const SpatialCell * cell, const bool calculateNonthermal, cuint popID, Real & rho ) {
Expand Down Expand Up @@ -1023,9 +1300,9 @@ namespace DRO {
return;
}

/*********
End velocity moment / thermal/non-thermal helper functions
*********/
/*********
End velocity moment / thermal/non-thermal helper functions
*********/

// Rho nonthermal:
VariableRhoNonthermal::VariableRhoNonthermal(cuint _popID): DataReductionOperator(),popID(_popID) {
Expand Down
Loading
Loading