Commit 5e6e7402 authored by Philipp Schäfer's avatar Philipp Schäfer
Browse files

Atmospheric Ray Tracing

- renamed CAdaptiveIntegrator to CAdaptiveSolver and moved it to own file (hidden outside project)
parent 656a91fc
......@@ -39,7 +39,6 @@ namespace ITAPropagationPathSim
{
private:
class CWorker;
class CAdaptiveIntegrator;
public:
SimulationSettings settings;
......
#include "AdaptiveSolver.h"
// ITA includes
#include <ITAException.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/ODESolver/ODESolver.h>
// STD
#include <cmath>
using namespace ITAPropagationPathSim::AtmosphericRayTracing;
void Simulation::CAdaptiveSolver::Process(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
if (rSettings.AdaptiveIntegration.bActive)
AdaptiveIntegrationStep(r, s, atmosphere);
else
IntegrationStep(r, s, atmosphere);
}
void Simulation::CAdaptiveSolver::IntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
{
switch (rSettings.solverMethod)
{
case SolverMethod::EULER:
ODESolver::Euler(r, s, dCurrentDt, atmosphere);
break;
case SolverMethod::RUNGE_KUTTA:
ODESolver::RungeKutta(r, s, dCurrentDt, atmosphere);
break;
default:
ITA_EXCEPT_INVALID_PARAMETER("Unknown solver method.");
}
}
void Simulation::CAdaptiveSolver::AdaptiveIntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
if (bIncreaseDtInNextStep)
TryIncreaseTimeStep();
VistaVector3D rNew;
VistaVector3D sNew;
while (1)
{
rNew = VistaVector3D(r);
sNew = VistaVector3D(s);
IntegrationStep(rNew, sNew, atmosphere);
const double rz = abs(r[Vista::Z]);
const double rzNew = abs(rNew[Vista::Z]);
const VistaVector3D nNew = ODESolver::SlownessToNormal(sNew);
const double error = EstimateIntegrationError(rz, rzNew, nNew, atmosphere);
if (error <= rSettings.AdaptiveIntegration.dMaxError)
{
bIncreaseDtInNextStep = error < rSettings.AdaptiveIntegration.dUncriticalError;
break;
}
if (!TryDecreaseTimeStep())
break;
}
r = rNew;
s = sNew;
UpdateTimeFramePortion();
}
double Simulation::CAdaptiveSolver::EstimateIntegrationError(const double& rz1, const double& rz2, const VistaVector3D& n2, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
{
const VistaVector3D windVectorError = atmosphere.WindVector(rz2) - atmosphere.WindVector(rz1) - atmosphere.WindVectorGradient(rz1) * (rz2 - rz1);
const double error = windVectorError[Vista::X] * n2[Vista::X] + windVectorError[Vista::Y] * n2[Vista::Y]; //Same as .Dot since z-component is zero
return std::abs(error);
}
bool Simulation::CAdaptiveSolver::TryIncreaseTimeStep()
{
if (iDtAdaptationLevel <= 0 || !IncreasedDtFitsInRestOfTimeFrame())
return false;
dCurrentDt *= 2;
iDtAdaptationLevel--;
return true;
}
bool Simulation::CAdaptiveSolver::TryDecreaseTimeStep()
{
if (iDtAdaptationLevel >= rSettings.AdaptiveIntegration.iMaxAdaptationLevel)
return false;
dCurrentDt /= 2;
iDtAdaptationLevel++;
return true;
}
void Simulation::CAdaptiveSolver::UpdateTimeFramePortion()
{
iTimeFramePortion += 1 << (rSettings.AdaptiveIntegration.iMaxAdaptationLevel - iDtAdaptationLevel);
if (iTimeFramePortion == maxDtPortion)
iTimeFramePortion = 0;
}
bool Simulation::CAdaptiveSolver::IncreasedDtFitsInRestOfTimeFrame() const
{
const unsigned int newDtPortion = 1 << (rSettings.AdaptiveIntegration.iMaxAdaptationLevel - (iDtAdaptationLevel - 1));
return ((maxDtPortion - iTimeFramePortion) % newDtPortion) == 0;
}
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_ADAPTIVESOLVER
#define IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_ADAPTIVESOLVER
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/Settings.h>
// ITA includes
#include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace Simulation
{
class CAdaptiveSolver
{
private:
const SimulationSettings rSettings;
double dCurrentDt = rSettings.dIntegrationTimeStep;
unsigned int iDtAdaptationLevel = 0;
const unsigned int maxDtPortion = 1 << rSettings.AdaptiveIntegration.iMaxAdaptationLevel; //! = 2 ^ iMaxAdaptionLevel
unsigned int iTimeFramePortion = 0; //! Portion of processed time of an unadapted time step. Integer between 0 and maxDtPortion.
bool bIncreaseDtInNextStep = false; //! Indicates whether the error is small enough that the timestep can be increased during next step.
public:
CAdaptiveSolver(const SimulationSettings& settings) : rSettings(settings) {}
public:
double GetCurrentStepSize() const { return dCurrentDt; }
void Process(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere);
private:
void IntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere) const;
void AdaptiveIntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere);
bool TryIncreaseTimeStep();
bool TryDecreaseTimeStep();
double EstimateIntegrationError(const double& rz1, const double& rz2, const VistaVector3D& n2, const ITAGeo::CStratifiedAtmosphere& atmosphere) const;
void UpdateTimeFramePortion();
bool IncreasedDtFitsInRestOfTimeFrame() const;
};
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_ADAPTIVESOLVER
\ No newline at end of file
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/Engine.h>
// ITA includes
#include <ITAException.h>
#include "AdaptiveSolver.h"
#include <ITAPropagationPathSim/AtmosphericRayTracing/ODESolver/ODESolver.h>
// STD
......@@ -9,124 +9,15 @@
using namespace ITAPropagationPathSim::AtmosphericRayTracing::Simulation;
class CEngine::CAdaptiveIntegrator
{
private:
const SimulationSettings rSettings;
double dCurrentDt = rSettings.dIntegrationTimeStep;
unsigned int iDtAdaptationLevel = 0;
const unsigned int maxDtPortion = 1 << rSettings.AdaptiveIntegration.iMaxAdaptationLevel; //! = 2 ^ iMaxAdaptionLevel
unsigned int iTimeFramePortion = 0; //! Portion of processed time of an unadapted time step. Integer between 0 and maxDtPortion.
bool bIncreaseDtInNextStep = false; //! Indicates whether the error is small enough that the timestep can be increased during next step.
public:
CAdaptiveIntegrator(const SimulationSettings& settings) : rSettings(settings) {}
private:
void IntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
{
switch (rSettings.solverMethod)
{
case SolverMethod::EULER:
ODESolver::Euler(r, s, dCurrentDt, atmosphere);
break;
case SolverMethod::RUNGE_KUTTA:
ODESolver::RungeKutta(r, s, dCurrentDt, atmosphere);
break;
default:
ITA_EXCEPT_INVALID_PARAMETER("Unknown solver method.");
}
}
void UpdateTimeFramePortion()
{
iTimeFramePortion += 1 << (rSettings.AdaptiveIntegration.iMaxAdaptationLevel - iDtAdaptationLevel);
if (iTimeFramePortion == maxDtPortion)
iTimeFramePortion = 0;
}
bool IncreasedDtFitsInRestOfTimeFrame() const
{
const unsigned int newDtPortion = 1 << (rSettings.AdaptiveIntegration.iMaxAdaptationLevel - (iDtAdaptationLevel - 1));
return ( (maxDtPortion-iTimeFramePortion) % newDtPortion) == 0;
}
bool TryIncreaseTimeStep()
{
if (iDtAdaptationLevel <= 0 || !IncreasedDtFitsInRestOfTimeFrame())
return false;
dCurrentDt *= 2;
iDtAdaptationLevel--;
return true;
}
bool TryDecreaseTimeStep()
{
if (iDtAdaptationLevel >= rSettings.AdaptiveIntegration.iMaxAdaptationLevel)
return false;
dCurrentDt /= 2;
iDtAdaptationLevel++;
return true;
}
double EstimateIntegrationError(const double& rz1, const double& rz2, const VistaVector3D& n2, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
{
const VistaVector3D windVectorError = atmosphere.WindVector(rz2) - atmosphere.WindVector(rz1) - atmosphere.WindVectorGradient(rz1) * (rz2 - rz1);
const double error = windVectorError[Vista::X] * n2[Vista::X] + windVectorError[Vista::Y] * n2[Vista::Y]; //Same as .Dot since z-component is zero
return std::abs(error);
}
void AdaptiveIntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
if (bIncreaseDtInNextStep)
TryIncreaseTimeStep();
VistaVector3D rNew;
VistaVector3D sNew;
while (1)
{
rNew = VistaVector3D(r);
sNew = VistaVector3D(s);
IntegrationStep(rNew, sNew, atmosphere);
const double rz = abs(r[Vista::Z]);
const double rzNew = abs(rNew[Vista::Z]);
const VistaVector3D nNew = ODESolver::SlownessToNormal(sNew);
const double error = EstimateIntegrationError(rz, rzNew, nNew, atmosphere);
if (error <= rSettings.AdaptiveIntegration.dMaxError)
{
bIncreaseDtInNextStep = error < rSettings.AdaptiveIntegration.dUncriticalError;
break;
}
if (!TryDecreaseTimeStep())
break;
}
r = rNew;
s = sNew;
UpdateTimeFramePortion();
}
public:
double GetCurrentStepSize() const { return dCurrentDt; }
void Process(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
if (rSettings.AdaptiveIntegration.bActive)
AdaptiveIntegrationStep(r, s, atmosphere);
else
IntegrationStep(r, s, atmosphere);
}
};
class CEngine::CWorker
{
private:
CRay& rRay;
const IAbortCriterion& rAbortCriterion; //!< Reference to externally defined abort criterion.
CAdaptiveIntegrator mIntegrator;
CAdaptiveSolver mSolver;
public:
CWorker(CRay& ray, const SimulationSettings& simSettings, const IAbortCriterion& abortCriterion) :
rRay(ray), rAbortCriterion(abortCriterion), mIntegrator(CAdaptiveIntegrator(simSettings)) {}
rRay(ray), rAbortCriterion(abortCriterion), mSolver(CAdaptiveSolver(simSettings)) {}
private:
bool ReflectionOccured(const double& rz1, const double& rz2) const
......@@ -154,9 +45,9 @@ class CEngine::CWorker
while ( !rAbortCriterion.AbortRequested(rRay) )
{
mIntegrator.Process(r, s, atmosphere);
mSolver.Process(r, s, atmosphere);
n = ODESolver::SlownessToNormal(s);
const double dt = mIntegrator.GetCurrentStepSize();
const double dt = mSolver.GetCurrentStepSize();
//if (ReflectionOccured(1, 1))
//{
......
......@@ -5,6 +5,8 @@ set( RelativeDir "src/ITAPropagationPathSim/AtmosphericRayTracing/Simulation" )
set( RelativeSourceGroup "Source Files\\ITAPropagationPathSim\\AtmosphericRayTracing\\Simulation" )
set( DirFiles
AdaptiveSolver.h
AdaptiveSolver.cpp
Engine.cpp
)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment