Skip to content

Oscillation suppression for DG #131

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
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Oscillation suppression for DG #131

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
2618b3c
Add infrastructure to export troubled cells
jbreue16 May 31, 2023
00c0d12
WIP oscillation suppression
jbreue16 Jan 28, 2024
55e0480
WIP add subcell limiting
jbreue16 Dec 8, 2023
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
34 changes: 34 additions & 0 deletions include/cadet/SolutionExporter.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,13 @@ class CADET_API ISolutionExporter
* @return @c true if particles are always represented by a single element, otherwise @c false
*/
virtual bool isParticleLumped() const CADET_NOEXCEPT = 0;

/**
* @brief Returns whether the associated model discretization has a smoothness indicator
* @details discretizations may contain a smoothness indicator for oscillation suppression
* @return @c true if indicator is present, otherwise @c false
*/
virtual bool hasSmoothnessIndicator() const CADET_NOEXCEPT = 0;

/**
* @brief Returns whether the primary coordinate is always a single element
Expand All @@ -85,13 +92,27 @@ class CADET_API ISolutionExporter
*/
virtual unsigned int numPrimaryCoordinates() const CADET_NOEXCEPT = 0;

/**
* @brief If existent (i.e. DG), returns the polynomial degree of the primary coordinate discretization.
* @details For axial flow columns discretized by DG, this is the axial polynomial degree
* @return Polynomial degree of primary coordinate discretization
*/
virtual unsigned int primaryPolynomialDegree() const CADET_NOEXCEPT = 0;

/**
* @brief Returns the number of secondary coordinates / points
* @details For axial flow columns, this is the number of radial points
* @return Number of secondary coordinates / points
*/
virtual unsigned int numSecondaryCoordinates() const CADET_NOEXCEPT = 0;

/**
* @brief If existent (i.e. DG), returns the polynomial degree of the secondary coordinate discretization.
* @details For axial flow columns discretized by DG, this is the radial polynomial degree
* @return Polynomial degree of secondary coordinate discretization
*/
virtual unsigned int secondaryPolynomialDegree() const CADET_NOEXCEPT = 0;

/**
* @brief Returns the number of inlet ports
* @return Number of inlet ports
Expand All @@ -117,6 +138,12 @@ class CADET_API ISolutionExporter
*/
virtual unsigned int numParticleShells(unsigned int parType) const CADET_NOEXCEPT = 0;

/**
* @brief If existent (i.e. DG), returns the polynomial degree of the particle discretization of a given particle type.
* @return Polynomial degree of particle discretization
*/
virtual unsigned int particlePolynomialDegree(unsigned int parType) const CADET_NOEXCEPT = 0;

/**
* @brief Returns the total number of bound states for all components of a given particle type
* @param [in] parType Particle type index
Expand Down Expand Up @@ -321,6 +348,13 @@ class CADET_API ISolutionExporter
*/
virtual int writeOutlet(double* buffer) const = 0;

/**
* @brief Writes the smoothness indicator if present, i.e. only for DG
* @details Writes the smoothness indicator for each DG element.
* @param [out] buffer Pointer to buffer that receives the data
* @return Number of written items
*/
virtual int writeSmoothnessIndicator(double* buffer) const = 0;

/**
* @brief Returns primary coordinates (e.g., axial for axial flow columns)
Expand Down
5 changes: 5 additions & 0 deletions include/common/Driver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,11 @@ void readDataOutputConfig(ParamProvider_t& pp, StorageConfig_t& cfg, const std::
cfg.storeVolume = pp.getBool("WRITE_" + dataType + "_VOLUME");
else
cfg.storeVolume = false;

if (pp.exists("WRITE_" + dataType + "_SMOOTHNESS_INDICATOR"))
cfg.storeSmoothnessIndicator = pp.getBool("WRITE_" + dataType + "_SMOOTHNESS_INDICATOR");
else
cfg.storeSmoothnessIndicator = false;
}

template <class ParamProvider_t>
Expand Down
105 changes: 71 additions & 34 deletions include/common/SolutionRecorderImpl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,14 +51,15 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
bool storeOutlet;
bool storeInlet;
bool storeVolume;
bool storeSmoothnessIndicator;
};

InternalStorageUnitOpRecorder() : InternalStorageUnitOpRecorder(UnitOpIndep) { }

InternalStorageUnitOpRecorder(UnitOpIdx idx) : _cfgSolution({false, false, false, true, false, false, false}),
_cfgSolutionDot({false, false, false, false, false, false, false}), _cfgSensitivity({false, false, false, true, false, false, false}),
_cfgSensitivityDot({false, false, false, true, false, false, false}), _storeTime(false), _storeCoordinates(false), _splitComponents(true), _splitPorts(true),
_singleAsMultiPortUnitOps(false), _keepBulkSingletonDim(true), _keepParticleSingletonDim(true), _curCfg(nullptr), _nComp(0), _nVolumeDof(0), _nAxialCells(0), _nRadialCells(0),
InternalStorageUnitOpRecorder(UnitOpIdx idx) : _cfgSolution({ false, false, false, true, false, false, false, false }),
_cfgSolutionDot({ false, false, false, false, false, false, false, false }), _cfgSensitivity({ false, false, false, true, false, false, false, false }),
_cfgSensitivityDot({ false, false, false, true, false, false, false, false }), _storeTime(false), _storeCoordinates(false),
_splitComponents(true), _splitPorts(true), _singleAsMultiPortUnitOps(false), _keepBulkSingletonDim(true), _keepParticleSingletonDim(true), _curCfg(nullptr), _nComp(0), _nVolumeDof(0), _nAxialPoints(0), _nRadialPoints(0),
_nInletPorts(0), _nOutletPorts(0), _numTimesteps(0), _numSens(0), _unitOp(idx), _needsReAlloc(false), _axialCoords(0), _radialCoords(0), _particleCoords(0)
{
}
Expand Down Expand Up @@ -127,8 +128,9 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
_keepBulkSingletonDim = exporter.hasPrimaryExtent();
_keepParticleSingletonDim = !exporter.isParticleLumped();

_nAxialCells = exporter.numPrimaryCoordinates();
_nRadialCells = exporter.numSecondaryCoordinates();
_nAxialPoints = exporter.numPrimaryCoordinates();
_nRadialPoints = exporter.numSecondaryCoordinates();
_axialPolyDeg = exporter.primaryPolynomialDegree();

// Query particle type specific structure
const unsigned int numParTypes = exporter.numParticleTypes();
Expand Down Expand Up @@ -158,7 +160,7 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
offset += _nParShells[i];
}
}

// Validate config
validateConfig(exporter, _cfgSolution);
validateConfig(exporter, _cfgSolutionDot);
Expand Down Expand Up @@ -246,6 +248,15 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
exporter.writeMobilePhase(v.data() + v.size() - sliceSize);
}

if (_curCfg->storeSmoothnessIndicator)
{
const int sliceSize = exporter.numMobilePhaseDofs() / (exporter.primaryPolynomialDegree() + 1u); // = nComp * nCells
std::vector<double>& v = _curStorage->smoothnessIndicator;

v.resize(v.size() + sliceSize);
exporter.writeSmoothnessIndicator(v.data() + v.size() - sliceSize);
}

if (_curCfg->storeParticle)
{
for (std::size_t parType = 0; parType < _nParShells.size(); ++parType)
Expand Down Expand Up @@ -475,6 +486,7 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
inline double const* solid(unsigned int parType = 0) const CADET_NOEXCEPT { return _data.solid[parType].data(); }
inline double const* flux() const CADET_NOEXCEPT { return _data.flux.data(); }
inline double const* volume() const CADET_NOEXCEPT { return _data.volume.data(); }
inline double const* smoothnessIndicator() const CADET_NOEXCEPT { return _data.smoothnessIndicator.data(); }
inline double const* inletDot() const CADET_NOEXCEPT { return _dataDot.inlet.data(); }
inline double const* outletDot() const CADET_NOEXCEPT { return _dataDot.outlet.data(); }
inline double const* bulkDot() const CADET_NOEXCEPT { return _dataDot.bulk.data(); }
Expand Down Expand Up @@ -510,6 +522,7 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
std::vector<std::vector<double>> solid;
std::vector<double> flux;
std::vector<double> volume;
std::vector<double> smoothnessIndicator;
};

inline void beginSensitivity(unsigned int sensIdx)
Expand All @@ -534,6 +547,7 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
cfg.storeSolid = exporter.hasSolidPhase() && cfg.storeSolid;
cfg.storeFlux = exporter.hasParticleFlux() && cfg.storeFlux;
cfg.storeVolume = exporter.hasVolume() && cfg.storeVolume;
cfg.storeSmoothnessIndicator = exporter.hasSmoothnessIndicator() && cfg.storeSmoothnessIndicator;
}

inline void allocateMemory(const ISolutionExporter& exporter)
Expand Down Expand Up @@ -568,6 +582,9 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder

if (_curCfg->storeVolume)
_curStorage->volume.reserve(nAllocTimesteps * exporter.numVolumeDofs());

if (_curCfg->storeSmoothnessIndicator)
_curStorage->smoothnessIndicator.reserve(nAllocTimesteps * _nComp * _nAxialPoints / (_axialPolyDeg + 1));
}

template <typename Writer_t>
Expand All @@ -590,8 +607,8 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
}
else
{
oss << prefix << "_OUTLET_PORT_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << port
<< "_COMP_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << comp;
oss << prefix << "_OUTLET_PORT_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << port
<< "_COMP_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << comp;
}

writer.template vector<double>(oss.str(), _numTimesteps, _curStorage->outlet.data() + comp + port * _nComp, _nComp * _nOutletPorts);
Expand Down Expand Up @@ -632,13 +649,13 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
oss << prefix << "_OUTLET";
if ((_nOutletPorts == 1) && !_singleAsMultiPortUnitOps)
{
const std::vector<std::size_t> layout = {_numTimesteps, _nComp};
const std::vector<std::size_t> layout = { _numTimesteps, _nComp };
debugCheckTensorLayout(layout, _curStorage->outlet.size());
writer.template tensor<double>(oss.str(), layout.size(), layout.data(), _curStorage->outlet.data());
}
else
{
const std::vector<std::size_t> layout = {_numTimesteps, _nOutletPorts, _nComp};
const std::vector<std::size_t> layout = { _numTimesteps, _nOutletPorts, _nComp };
debugCheckTensorLayout(layout, _curStorage->outlet.size());
writer.template tensor<double>(oss.str(), layout.size(), layout.data(), _curStorage->outlet.data());
}
Expand All @@ -663,8 +680,8 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
}
else
{
oss << prefix << "_INLET_PORT_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << port
<< "_COMP_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << comp;
oss << prefix << "_INLET_PORT_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << port
<< "_COMP_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << comp;
}

writer.template vector<double>(oss.str(), _numTimesteps, _curStorage->inlet.data() + comp + port * _nComp, _nComp * _nInletPorts);
Expand Down Expand Up @@ -705,13 +722,13 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
oss << prefix << "_INLET";
if ((_nInletPorts == 1) && !_singleAsMultiPortUnitOps)
{
const std::vector<std::size_t> layout = {_numTimesteps, _nComp};
const std::vector<std::size_t> layout = { _numTimesteps, _nComp };
debugCheckTensorLayout(layout, _curStorage->inlet.size());
writer.template tensor<double>(oss.str(), layout.size(), layout.data(), _curStorage->inlet.data());
}
else
{
const std::vector<std::size_t> layout = {_numTimesteps, _nInletPorts, _nComp};
const std::vector<std::size_t> layout = { _numTimesteps, _nInletPorts, _nComp };
debugCheckTensorLayout(layout, _curStorage->inlet.size());
writer.template tensor<double>(oss.str(), layout.size(), layout.data(), _curStorage->inlet.data());
}
Expand All @@ -728,10 +745,10 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
layout.reserve(4);
layout.push_back(_numTimesteps);

if ((_keepBulkSingletonDim && (_nAxialCells == 1)) || (_nAxialCells > 1))
layout.push_back(_nAxialCells);
if (_nRadialCells > 0)
layout.push_back(_nRadialCells);
if ((_keepBulkSingletonDim && (_nAxialPoints == 1)) || (_nAxialPoints > 1))
layout.push_back(_nAxialPoints);
if (_nRadialPoints > 0)
layout.push_back(_nRadialPoints);
layout.push_back(_nComp);

debugCheckTensorLayout(layout, _curStorage->bulk.size());
Expand All @@ -745,10 +762,10 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
layout.reserve(5);
layout.push_back(_numTimesteps);

if ((_keepBulkSingletonDim && (_nAxialCells == 1)) || (_nAxialCells > 1))
layout.push_back(_nAxialCells);
if (_nRadialCells > 0)
layout.push_back(_nRadialCells);
if ((_keepBulkSingletonDim && (_nAxialPoints == 1)) || (_nAxialPoints > 1))
layout.push_back(_nAxialPoints);
if (_nRadialPoints > 0)
layout.push_back(_nRadialPoints);

if (_nParShells.size() <= 1)
{
Expand Down Expand Up @@ -805,10 +822,10 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
layout.reserve(5);
layout.push_back(_numTimesteps);

if ((_keepBulkSingletonDim && (_nAxialCells == 1)) || (_nAxialCells > 1))
layout.push_back(_nAxialCells);
if (_nRadialCells > 0)
layout.push_back(_nRadialCells);
if ((_keepBulkSingletonDim && (_nAxialPoints == 1)) || (_nAxialPoints > 1))
layout.push_back(_nAxialPoints);
if (_nRadialPoints > 0)
layout.push_back(_nRadialPoints);

if (_nParShells.size() <= 1)
{
Expand Down Expand Up @@ -865,10 +882,10 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder

layout.push_back(_numTimesteps);
layout.push_back(_nParShells.size());
if ((_keepBulkSingletonDim && (_nAxialCells == 1)) || (_nAxialCells > 1))
layout.push_back(_nAxialCells);
if (_nRadialCells > 0)
layout.push_back(_nRadialCells);
if ((_keepBulkSingletonDim && (_nAxialPoints == 1)) || (_nAxialPoints > 1))
layout.push_back(_nAxialPoints);
if (_nRadialPoints > 0)
layout.push_back(_nRadialPoints);
layout.push_back(_nComp);

debugCheckTensorLayout(layout, _curStorage->flux.size());
Expand All @@ -884,6 +901,24 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder
oss << prefix << "_VOLUME";
writer.template matrix<double>(oss.str(), _numTimesteps, _nVolumeDof, _curStorage->volume.data(), 1);
}

if (_curCfg->storeSmoothnessIndicator)
{
oss.str("");
oss << prefix << "_SMOOTHNESS_INDICATOR";

std::vector<std::size_t> layout(0);
layout.reserve(4);
layout.push_back(_numTimesteps);

if (_nAxialPoints > 0)
layout.push_back(_nAxialPoints / (_axialPolyDeg + 1));
layout.push_back(_nComp);

debugCheckTensorLayout(layout, _curStorage->smoothnessIndicator.size());

writer.template tensor<double>(oss.str(), layout.size(), layout.data(), _curStorage->smoothnessIndicator.data());
}
}

inline void clear(Storage& s)
Expand All @@ -900,6 +935,7 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder

s.flux.clear();
s.volume.clear();
s.smoothnessIndicator.clear();
}

template <typename T>
Expand Down Expand Up @@ -945,8 +981,9 @@ class InternalStorageUnitOpRecorder : public ISolutionRecorder

unsigned int _nComp;
unsigned int _nVolumeDof;
unsigned int _nAxialCells;
unsigned int _nRadialCells;
unsigned int _nAxialPoints;
unsigned int _nRadialPoints;
int _axialPolyDeg;
unsigned int _nInletPorts;
unsigned int _nOutletPorts;
std::vector<unsigned int> _nParShells;
Expand Down Expand Up @@ -1119,7 +1156,7 @@ class InternalStorageSystemRecorder : public ISolutionRecorder
}
}
}

template <typename Writer_t>
void writeSolution(Writer_t& writer)
{
Expand Down
Loading
Loading