Skip to content

Commit

Permalink
Flow deltaT adjustment on thin volumes (#3045)
Browse files Browse the repository at this point in the history 
* fix #3008

* format

* minor

* address review comments
  • Loading branch information
@shaomeng
shaomeng committed Feb 23, 2022
1 parent 809d4d6 commit 524915a
Show file tree
Hide file tree
Showing 6 changed files with 61 additions and 20 deletions.
3 changes: 3 additions & 0 deletions include/vapor/Advection.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,9 @@ class FLOW_API Advection final {
// Returns the value after adjustment.
float _applyPeriodic(float val, float min, float max) const;

// Print return code if it's non-zero and compiled in debug mode.
void _printNonZero(int rtn, const char *file, const char *func, int line) const;

void _calculateParticleIntegratedValue(Particle &p, const Particle &prev, const Field *scalarField, const bool skipNonZero, const float distScale,
const std::vector<double> &integrateWithinVolumeMin, const std::vector<double> &integrateWithinVolumeMax) const;
static bool _isParticleInsideVolume(const Particle &p, const std::vector<double> &min, const std::vector<double> &max);
Expand Down
2 changes: 1 addition & 1 deletion include/vapor/FlowRenderer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ class RENDER_API FlowRenderer final : public Renderer {
double newTime) const; // New time to assign to particles

// Print return code if it's non-zero and compiled in debug mode.
void _printNonZero(int rtn, const char *file, const char *func, int line);
void _printNonZero(int rtn, const char *file, const char *func, int line) const;

}; // End of class FlowRenderer

Expand Down
43 changes: 34 additions & 9 deletions lib/flow/Advection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,14 @@ int Advection::AdvectSteps(Field *velocity, double deltaT, size_t maxSteps, ADVE
Particle p1;
int rv = 0;
switch (method) {
case ADVECTION_METHOD::EULER: rv = _advectEuler(velocity, past0, dt, p1); break;
case ADVECTION_METHOD::RK4: rv = _advectRK4(velocity, past0, dt, p1); break;
case ADVECTION_METHOD::EULER:
rv = _advectEuler(velocity, past0, dt, p1);
_printNonZero(rv, __FILE__, __func__, __LINE__);
break;
case ADVECTION_METHOD::RK4:
rv = _advectRK4(velocity, past0, dt, p1);
_printNonZero(rv, __FILE__, __func__, __LINE__);
break;
}

if (rv == 0) { // Advection successful!
Expand Down Expand Up @@ -239,8 +245,14 @@ int Advection::AdvectTillTime(Field *velocity, double startT, double deltaT, dou
Particle p1;
int rv = 0;
switch (method) {
case ADVECTION_METHOD::EULER: rv = _advectEuler(velocity, p0, dt, p1); break;
case ADVECTION_METHOD::RK4: rv = _advectRK4(velocity, p0, dt, p1); break;
case ADVECTION_METHOD::EULER:
rv = _advectEuler(velocity, p0, dt, p1);
_printNonZero(rv, __FILE__, __func__, __LINE__);
break;
case ADVECTION_METHOD::RK4:
rv = _advectRK4(velocity, p0, dt, p1);
_printNonZero(rv, __FILE__, __func__, __LINE__);
break;
}
if (rv != 0) { // Advection wasn't successful for some reason...
break;
Expand Down Expand Up @@ -519,18 +531,22 @@ int Advection::_advectEuler(Field *velocity, const Particle &p0, double dt, Part

int Advection::_advectRK4(Field *velocity, const Particle &p0, double dt, Particle &p1) const
{
glm::vec3 k1, k2, k3, k4;
double dt_half = dt * 0.5;
float dt32 = float(dt); // glm is strict about data types (which is a good thing).
float dt_half32 = float(dt_half); // glm is strict about data types (which is a good thing).
int rv;
glm::vec3 k1, k2, k3, k4;
const double dt_half = dt * 0.5;
const float dt32 = float(dt); // glm is strict about data types (which is a good thing).
const float dt_half32 = float(dt_half); // glm is strict about data types (which is a good thing).
int rv = 0;
rv = velocity->GetVelocity(p0.time, p0.location, k1);
_printNonZero(rv, __FILE__, __func__, __LINE__);
if (rv != 0) return rv;
rv = velocity->GetVelocity(p0.time + dt_half, p0.location + dt_half32 * k1, k2);
_printNonZero(rv, __FILE__, __func__, __LINE__);
if (rv != 0) return rv;
rv = velocity->GetVelocity(p0.time + dt_half, p0.location + dt_half32 * k2, k3);
_printNonZero(rv, __FILE__, __func__, __LINE__);
if (rv != 0) return rv;
rv = velocity->GetVelocity(p0.time + dt, p0.location + dt32 * k3, k4);
_printNonZero(rv, __FILE__, __func__, __LINE__);
if (rv != 0) return rv;
p1.location = p0.location + dt32 / 6.0f * (k1 + 2.0f * (k2 + k3) + k4);
p1.time = p0.time + dt;
Expand Down Expand Up @@ -656,6 +672,15 @@ float Advection::_applyPeriodic(float val, float min, float max) const
}
}

void Advection::_printNonZero(int rtn, const char *file, const char *func, int line) const
{
#ifndef NDEBUG
if (rtn != 0) { // only print non-zero values
printf("Rtn == %d: %s:(%s):%d\n", rtn, file, func, line);
}
#endif
}

auto Advection::GetValueVarName() const -> std::string { return _valueVarName; }

auto Advection::GetPropertyVarNames() const -> std::vector<std::string> { return _propertyVarNames; }
Expand Down
4 changes: 1 addition & 3 deletions lib/flow/Particle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,9 +56,7 @@ void Particle::SetSpecial(bool isSpecial)
time = 0.0f;
value = 0.0f;
}
location.x = 0.0f;
location.y = 0.0f;
location.z = 0.0f;
location = {0.0f, 0.0f, 0.0f};
}

bool Particle::IsSpecial() const { return (std::isnan(time) && std::isnan(value)); }
21 changes: 17 additions & 4 deletions lib/flow/VaporField.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,13 +238,15 @@ int VaporField::GetVelocity(double time, const glm::vec3 &pos, glm::vec3 &veloci
}
auto hasMissing = glm::equal(velocity, missingV);
float mult = _params_locked ? _c_vel_mult : _params->GetVelocityMultiplier();
if (glm::any(hasMissing))
if (glm::any(hasMissing)) {
return MISSING_VAL;
} //
else {
velocity *= mult;
return 0;
return SUCCESS;
}
} else {
} // Finish steady case
else {
float mult = _params->GetVelocityMultiplier();

// First check if the query time is within range
Expand Down Expand Up @@ -441,7 +443,7 @@ int VaporField::CalcDeltaTFromCurrentTimeStep(double &delT) const
return rv;
else {
auto mag = glm::length(vel);
if (mag > maxmag) maxmag = mag;
maxmag = std::max(maxmag, mag);
}
}

Expand Down Expand Up @@ -469,6 +471,17 @@ int VaporField::CalcDeltaTFromCurrentTimeStep(double &delT) const
const double actualNum = double(glm::distance(minxyz, maxxyz)) / double(maxmag);
delT = actualNum / desiredNum;

// Another logic in determing the size of deltaT:
// if deltaT will send a particle out of the domain in just one step in any direction,
// then we halve deltaT until the particle can go one step inside of the volume.
float smallestD = std::numeric_limits<float>::max();
if (!VelocityNames[0].empty()) smallestD = std::min(smallestD, std::abs(minxyz.x - maxxyz.x));
if (!VelocityNames[1].empty()) smallestD = std::min(smallestD, std::abs(minxyz.y - maxxyz.y));
if (!VelocityNames[2].empty()) smallestD = std::min(smallestD, std::abs(minxyz.z - maxxyz.z));
while (maxmag * delT > double(smallestD)) { //
delT /= 2.0;
}

return 0;
}

Expand Down
8 changes: 5 additions & 3 deletions lib/render/FlowRenderer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,14 +278,16 @@ int FlowRenderer::_paintGL(bool fast)

Progress::StartIndefinite("Performing flowline calculations");
Progress::Update(0);
_advection.AdvectSteps(&_velocityField, deltaT, numOfSteps);
rv = _advection.AdvectSteps(&_velocityField, deltaT, numOfSteps);
_printNonZero(rv, __FILE__, __func__, __LINE__);

// If the advection is bi-directional
if (_2ndAdvection) {
assert(deltaT > 0.0);
auto deltaT2 = deltaT * -1.0;

_2ndAdvection->AdvectSteps(&_velocityField, deltaT2, numOfSteps);
rv = _2ndAdvection->AdvectSteps(&_velocityField, deltaT2, numOfSteps);
_printNonZero(rv, __FILE__, __func__, __LINE__);
}
Progress::Finish();
}
Expand Down Expand Up @@ -1280,7 +1282,7 @@ int FlowRenderer::_updateAdvectionPeriodicity(flow::Advection *advc)
return 0;
}

void FlowRenderer::_printNonZero(int rtn, const char *file, const char *func, int line)
void FlowRenderer::_printNonZero(int rtn, const char *file, const char *func, int line) const
{
#ifndef NDEBUG
if (rtn != 0) { // only print non-zero values
Expand Down

0 comments on commit 524915a

Please sign in to comment.