i use the library to write a termNormal but it the data use matlab to show the 3d object it not correct
Closed this issue · 1 comments
YuDe95 commented
the code as follows,although it not very clean
- main.cpp
// burnNormal.cpp : containt "main" function。
//
#define _USE_MATH_DEFINES
#include <levelset/levelset.hpp>
#include <cmath>
#include <numeric>
#include <functional>
#include <cfloat>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <cstring>
#include "BurnNormalSchemeData.h"
#include <opencv2/opencv.hpp>
#define VISUALIZE_BY_OPENCV
#include <helperOC/helperOC.hpp>
#include <helperOC/DynSys/DynSys/DynSysSchemeData.hpp>
#include <helperOC/DynSys/Quad4D/Quad4D.hpp>
#include <helperOC/ValFuncs/proj.hpp>
#include <helperOC/ValFuncs/visSetIm.hpp>
typedef enum ApproximationAccuracy_Type {
ApproximationAccuracy_Invalid,
ApproximationAccuracy_low,
ApproximationAccuracy_medium,
ApproximationAccuracy_high,
ApproximationAccuracy_veryHigh,
}ApproximationAccuracy_Type;
#include <stdio.h>
#include <cuda_runtime.h>
#include <cuda.h>
#include <iostream>
int main(int argc, char* argv[])
{
//bool debug_dump_file = false;
bool debug_dump_file = true;
bool dump_file = false;
if (argc >= 2) {
dump_file = (atoi(argv[1]) == 0) ? false : true;
}
size_t line_length_of_chunk = 1;
if (argc >= 3) {
line_length_of_chunk = atoi(argv[2]);
}
bool useCuda = false;
if (argc >= 4) {
useCuda = (atoi(argv[3]) == 0) ? false : true;
}
int num_of_threads = 0;
if (argc >= 5) {
num_of_threads = atoi(argv[4]);
}
int num_of_gpus = 0;
if (argc >= 6) {
num_of_gpus = atoi(argv[5]);
}
levelset::DelayedDerivMinMax_Type delayedDerivMinMax = levelset::DelayedDerivMinMax_Disable;
if (argc >= 7) {
switch (atoi(argv[6])) {
default:
case 0:
delayedDerivMinMax = levelset::DelayedDerivMinMax_Disable;
break;
case 1:
delayedDerivMinMax = levelset::DelayedDerivMinMax_Always;
break;
case 2:
delayedDerivMinMax = levelset::DelayedDerivMinMax_Adaptive;
break;
}
}
bool enable_user_defined_dynamics_on_gpu = true;
if (argc >= 8) {
enable_user_defined_dynamics_on_gpu = (atoi(argv[7]) == 0) ? false : true;
}
const FLOAT_TYPE tMax = (FLOAT_TYPE)15; //!< End time.
// const FLOAT_TYPE tMax = 0.35; //!< End time.
// const FLOAT_TYPE tMax = 0.1; //!< End time.
const int plotSteps = 9; //!< How many intermediate plots to produce?
// const int plotSteps = 2; //!< How many intermediate plots to produce?
const FLOAT_TYPE t0 = 0.0; //!< Start time.
FLOAT_TYPE tPlot = (tMax - t0) / (plotSteps - 1);
//! How close (relative) do we need to get to tMax to be considered finished?
const FLOAT_TYPE small_ratio = 100.; //!<
const FLOAT_TYPE eps = std::numeric_limits<FLOAT_TYPE>::epsilon(); //!<
FLOAT_TYPE small = small_ratio * eps; //!<
//! Problem Parameters.
//! Radius of target circle(positive).
const FLOAT_TYPE targetRadius = 5;
//! Speed of the evader(positive constant).
size_t num_of_dimensions = 3;
//beacls::FloatVec mins{ -6,-10,0 };
//beacls::FloatVec maxs{ +20,+10,(FLOAT_TYPE)(+2 * M_PI) };
//size_t Nx = 51;
//beacls::IntegerVec Ns(num_of_dimensions);
//Ns = { Nx, (size_t)std::ceil(Nx * (maxs[1] - mins[1]) / (maxs[0] - mins[0])), Nx - 1 };
//maxs[2] = (FLOAT_TYPE)(maxs[2] * (1 - 1. / Ns[2]));
beacls::FloatVec mins{ -56,-50,-40 };
beacls::FloatVec maxs{ +40,+40,(FLOAT_TYPE)(+20 * M_PI) };
size_t Nx = 100;
beacls::IntegerVec Ns(num_of_dimensions);
Ns = { Nx, (size_t)std::ceil(Nx * (maxs[1] - mins[1]) / (maxs[0] - mins[0])), Nx - 1 };
maxs[2] = (FLOAT_TYPE)(maxs[2] * (1 - 1. / Ns[2]));
//beacls::FloatVec mins{ 158, - 62, - 24 };
//beacls::FloatVec maxs{ 679 ,90 ,133 };
//beacls::IntegerVec Ns(num_of_dimensions);
//Ns = { 53 ,152, 157 };
levelset::ShapeCylinder* shape = new levelset::ShapeCylinder(
beacls::IntegerVec{ 2 },
beacls::FloatVec{ 0.,0.,0 },
targetRadius);
levelset::AddGhostExtrapolate* addGhostExtrapolate = new levelset::AddGhostExtrapolate();
std::vector<levelset::BoundaryCondition*> boundaryConditions(3);
boundaryConditions[0] = addGhostExtrapolate;
boundaryConditions[1] = addGhostExtrapolate;
boundaryConditions[2] = addGhostExtrapolate;
levelset::HJI_Grid* hJI_Grid = new levelset::HJI_Grid(
num_of_dimensions);
hJI_Grid->set_mins(mins);
hJI_Grid->set_maxs(maxs);
hJI_Grid->set_boundaryConditions(boundaryConditions);
hJI_Grid->set_Ns(Ns);
if (!hJI_Grid->processGrid()) {
return -1;
}
const beacls::UVecType type = useCuda ? beacls::UVecType_Cuda : beacls::UVecType_Vector;
ApproximationAccuracy_Type accuracy = ApproximationAccuracy_veryHigh;
// ApproximationAccuracy_Type accuracy = ApproximationAccuracy_high;
// ApproximationAccuracy_Type accuracy = ApproximationAccuracy_veryHigh;
// Choose spatial derivative approimation at appropriate level of accuracy.
levelset::SpatialDerivative* spatialDerivative = NULL;
switch (accuracy) {
case ApproximationAccuracy_low:
spatialDerivative = new levelset::UpwindFirstFirst(hJI_Grid, type);
break;
case ApproximationAccuracy_medium:
spatialDerivative = new levelset::UpwindFirstENO2(hJI_Grid, type);
break;
case ApproximationAccuracy_high:
spatialDerivative = new levelset::UpwindFirstENO3(hJI_Grid, type);
break;
case ApproximationAccuracy_veryHigh:
spatialDerivative = new levelset::UpwindFirstWENO5(hJI_Grid, type);
break;
default:
printf("Unkown accuracy level %d\n", accuracy);
return -1;
}
std::vector<levelset::PostTimestep_Exec_Type*> postTimestep_Execs;
levelset::Integrator* integrator;
FLOAT_TYPE factor_cfl = (FLOAT_TYPE)0.75;
FLOAT_TYPE max_step = (FLOAT_TYPE)8.0e16;
bool single_step = false;
// bool single_step = true;
bool stats = true;
levelset::TerminalEvent_Exec_Type* terminalEvent_Exec_Type = NULL;
BurnNormalSchemeData* schemeData = new BurnNormalSchemeData();
schemeData->set_spatialDerivative(spatialDerivative);
schemeData->set_grid(hJI_Grid);
levelset::TermNormal* schemeFunc = new levelset::TermNormal(schemeData,type);
// Choose integration approimation at appropriate level of accuracy.
switch (accuracy) {
case ApproximationAccuracy_low:
integrator = new levelset::OdeCFL1(schemeFunc, factor_cfl, max_step, postTimestep_Execs, single_step, stats, terminalEvent_Exec_Type);
break;
case ApproximationAccuracy_medium:
integrator = new levelset::OdeCFL2(schemeFunc, factor_cfl, max_step, postTimestep_Execs, single_step, stats, terminalEvent_Exec_Type);
break;
case ApproximationAccuracy_high:
integrator = new levelset::OdeCFL3(schemeFunc, factor_cfl, max_step, postTimestep_Execs, single_step, stats, terminalEvent_Exec_Type);
break;
case ApproximationAccuracy_veryHigh:
integrator = new levelset::OdeCFL3(schemeFunc, factor_cfl, max_step, postTimestep_Execs, single_step, stats, terminalEvent_Exec_Type);
break;
default:
printf("Unkown accuracy level %d\n", accuracy);
return -1;
}
beacls::FloatVec data;
//std::ifstream file;
//file.open(R"(D:\LocalGitCode\SRM\Test\xxxx.sdf)");
//if (!file.is_open()) {
// //cerr << "open error!" << endl;
// exit(1);
//}
//std::string line;
//for (int skipLine = 0; skipLine < 3; ++skipLine)
//{
// std::getline(file, line);
//}
//while (!file.eof())
//{
// std::getline(file, line);
// if (!line.empty())
// {
// data.push_back(std::stof(line));
// }
// //data.push_back(std::stof(line));
//}
shape->execute(hJI_Grid, data);
if (false) {
beacls::MatFStream* fs = beacls::openMatFStream(std::string("initial.mat"), beacls::MatOpenMode_Write);
save_vector(data, std::string("y0"), Ns, false, fs);
beacls::closeMatFStream(fs);
// beacls::FloatVec reload_data;
// load_vector(std::string("initial.mat"), reload_data);
// HJI_Grid* tmp_grid = new levelset::HJI_Grid();
// tmp_grid->load_grid(std::string("air3D_0_g_v7_3.mat"));
// tmp_grid->save_grid(std::string("air3D_0_g_v7_3_new.mat"),std::string("grid"));
// delete tmp_grid;
}
FLOAT_TYPE tNow = t0;
beacls::FloatVec y;
beacls::FloatVec y0;
std::cout << "data size:" << data.size();
while ((tMax - tNow) > small * tMax) {
y0 = data;
beacls::FloatVec tspan(2);
tspan[0] = tNow;
tspan[1] = HjiMin(tMax, tNow + tPlot);
if (debug_dump_file) {
std::stringstream ss;
ss << std::setprecision(5) << tNow << std::resetiosflags(std::ios_base::floatfield);
std::string filename = "normal_" + ss.str() + ".txt";
dump_vector(filename.c_str(), data);
}
tNow = integrator->execute(
y, tspan, y0, schemeData,
line_length_of_chunk, num_of_threads, num_of_gpus,
delayedDerivMinMax, enable_user_defined_dynamics_on_gpu);
data = y;
printf("tNow = %f\n", tNow);
////!< Visualize
//levelset::HJI_Grid* g2Dp;
//beacls::FloatVec data2Dp;
//g2Dp = helperOC::proj(data2Dp, schemeData->get_grid(), data, beacls::IntegerVec{ 0,1,1 });
//std::string file = R"(D:\LocalGitCode\SRM\SRMbase1.21\beacls\beacls-master\sources\samples\air3D\x64\Release\)";
//cv::Mat image_mat;
//helperOC::visSetIm(image_mat, image_mat, schemeData->get_grid(), data2Dp, std::vector<float>{0.0,1.0,1.0}, beacls::FloatVec(),true, file+std::to_string(tNow)+ ".png",cv::Size());
////cv::namedWindow(__func__, 0);
////cv::imshow(__func__, image_mat);
////cv::waitKey(1);
//cv::imshow(file, image_mat);
//cv::waitKey(1);
}
if (debug_dump_file) {
std::stringstream ss;
ss << std::setprecision(5) << tNow << std::resetiosflags(std::ios_base::floatfield);
std::string filename = "normal_" + ss.str() + ".txt";
dump_vector(filename.c_str(), data);
}
if (schemeData) delete schemeData;
if (schemeFunc) delete schemeFunc;
if (integrator) delete integrator;
if (hJI_Grid) delete hJI_Grid;
if (shape) delete shape;
if (addGhostExtrapolate) delete addGhostExtrapolate;
if (spatialDerivative) delete spatialDerivative;
}
- BurnNormalSchemeData.hpp
#pragma once
#include <levelset/levelset.hpp>
#include <cstdint>
#include <vector>
#include <cstddef>
#include <utility>
#include <Core/UVec.hpp>
using namespace std::rel_ops;
enum VelocityType
{
Velocity_Float,//velocity type is a single float value
Velocity_UVec//velocity type is a vector value use it to set diffrent velocity in surface
};
class BurnNormalSchemeData : public levelset::SchemeData
{
private:
FLOAT_TYPE velocity = 10.0;
beacls::UVec velocity_uvec;
FLOAT_TYPE input;
VelocityType velocity_type;
public:
BurnNormalSchemeData() :SchemeData(), velocity(10.0), input(0.0)
{
velocity = 10.0;
}
~BurnNormalSchemeData(){}
bool operator==(const BurnNormalSchemeData& rhs) const;
bool operator==(const SchemeData& rhs)const;
BurnNormalSchemeData* clone()const
{
return new BurnNormalSchemeData(*this);
}
bool hamFunc(
beacls::UVec& hamValue_uvec,
const FLOAT_TYPE t,
const beacls::UVec& data,
const std::vector<beacls::UVec>& derivs,
const size_t begin_index,
const size_t length
)const override;
bool normal3DHamFunc(
beacls::UVec& hamValue_uvec,
beacls::UVec& magnitude_uvec,
beacls::FloatVec& step_bound_invs,
const std::vector<beacls::UVec>& derivs_l,
const std::vector<beacls::UVec>& derivs_r,
const size_t begin_index,
const size_t length
) const override;
void setVelocity(FLOAT_TYPE velocity);
private:
/**
* Disable operator=.
*/
BurnNormalSchemeData& operator=(const BurnNormalSchemeData& rhs);
BurnNormalSchemeData(const BurnNormalSchemeData& rhs):
SchemeData(rhs),
velocity(0),
input(0)
{
}
/**
* @brief Determine the upwind direction,Either both sides agree in sign (take direction in which they agree),
* or characteristics are converging (take larger magnitude direction)..
*
* \param prodL
* \param prodR
* \param magL
* \param magR
* \return
*/
int calaFlowLValue(FLOAT_TYPE prodL, FLOAT_TYPE prodR, FLOAT_TYPE magL, FLOAT_TYPE magR) const;
/**
* @brief Determine the upwind direction,Either both sides agree in sign (take direction in which they agree),
* or characteristics are converging (take larger magnitude direction).
*
* \param prodL
* \param prodR
* \param magL
* \param magR
* \return
*/
int calaFlowRValue(FLOAT_TYPE prodL, FLOAT_TYPE prodR, FLOAT_TYPE magL, FLOAT_TYPE magR) const;
};
BurnNormalSchemeData.cpp###
#include "BurnNormalSchemeData.h"
//#include <algorithm>
//#include <typeinfo>
//#include <macro.hpp>
#include <levelset/Grids/HJI_Grid.hpp>
bool BurnNormalSchemeData::operator==(const BurnNormalSchemeData& rhs) const
{
if (this == &rhs)
{
return true;
}
else if (!SchemeData::operator==(rhs))
{
return false;
}
else if (velocity != rhs.velocity)
{
return false;
}
else if (input != rhs.input)
{
return false;
}
else
{
return true;
}
}
bool BurnNormalSchemeData::operator==(const SchemeData& rhs) const
{
if (this == &rhs)
{
return true;
}
else if (typeid(*this) != typeid(rhs))
{
return false;
}
else
{
return operator==(dynamic_cast<const BurnNormalSchemeData&>(rhs));
}
}
bool BurnNormalSchemeData::hamFunc(
beacls::UVec& hamValue_uvec,
const FLOAT_TYPE t,
const beacls::UVec&,
const std::vector<beacls::UVec>& derivs,
const size_t begin_index,
const size_t length) const
{
return false;
}
bool BurnNormalSchemeData::normal3DHamFunc(
beacls::UVec& hamValue_uvec,
beacls::UVec& magnitude_uvec,
beacls::FloatVec& step_bound_invs,
const std::vector<beacls::UVec>& derivs_l,
const std::vector<beacls::UVec>& derivs_r,
const size_t begin_index, const size_t length) const
{
//printf("tNow = %u\n", derivs_l[0].size());
beacls::reallocateAsSrc(magnitude_uvec, derivs_l[0]);
FLOAT_TYPE* magnitudeValue = beacls::UVec_<FLOAT_TYPE>(magnitude_uvec).ptr();
beacls::reallocateAsSrc(hamValue_uvec, derivs_l[0]);
FLOAT_TYPE* hamValue = beacls::UVec_<FLOAT_TYPE>(hamValue_uvec).ptr();
const FLOAT_TYPE* deriv0_l = beacls::UVec_<FLOAT_TYPE>(derivs_l[0]).ptr();
const FLOAT_TYPE* deriv1_l = beacls::UVec_<FLOAT_TYPE>(derivs_l[1]).ptr();
const FLOAT_TYPE* deriv2_l = beacls::UVec_<FLOAT_TYPE>(derivs_l[2]).ptr();
const FLOAT_TYPE* deriv0_r = beacls::UVec_<FLOAT_TYPE>(derivs_r[0]).ptr();
const FLOAT_TYPE* deriv1_r = beacls::UVec_<FLOAT_TYPE>(derivs_r[1]).ptr();
const FLOAT_TYPE* deriv2_r = beacls::UVec_<FLOAT_TYPE>(derivs_r[2]).ptr();
FLOAT_TYPE prodL0;
FLOAT_TYPE prodL1;
FLOAT_TYPE prodL2;
FLOAT_TYPE prodR0;
FLOAT_TYPE prodR1;
FLOAT_TYPE prodR2;
FLOAT_TYPE magL0;
FLOAT_TYPE magL1;
FLOAT_TYPE magL2;
FLOAT_TYPE magR0;
FLOAT_TYPE magR1;
FLOAT_TYPE magR2;
int flowL0;
int flowL1;
int flowL2;
int flowR0;
int flowR1;
int flowR2;
const levelset::HJI_Grid* hji_grid = get_grid();
auto dimensions = hji_grid->get_num_of_dimensions();
if (step_bound_invs.size() != hji_grid->get_num_of_dimensions())
{
step_bound_invs.resize(hji_grid->get_num_of_dimensions());
}
FLOAT_TYPE stepBoundMaxInv0 = 0.0;
FLOAT_TYPE stepBoundMaxInv1 = 0.0;
FLOAT_TYPE stepBoundMaxInv2 = 0.0;
for (size_t i = 0; i < length; ++i)
{
/*for (int dimension = 0; dimension < dimensions; ++dimension)
{
}*/
FLOAT_TYPE velocitytep = -1.0;
FLOAT_TYPE deriv0_l_i = deriv0_l[i];
FLOAT_TYPE deriv1_l_i = deriv1_l[i];
FLOAT_TYPE deriv2_l_i = deriv2_l[i];
FLOAT_TYPE deriv0_r_i = deriv0_r[i];
FLOAT_TYPE deriv1_r_i = deriv1_r[i];
FLOAT_TYPE deriv2_r_i = deriv2_r[i];
/* prodL0 = velocity * deriv0_l_i; prodL1 = velocity * deriv1_l_i; prodL2 = velocity * deriv2_l_i;
prodR0 = velocity * deriv0_r_i; prodR1 = velocity * deriv1_r_i; prodR2 = velocity * deriv2_r_i;*/
prodL0 = velocitytep * deriv0_l_i; prodL1 = velocitytep * deriv1_l_i; prodL2 = velocitytep * deriv2_l_i;
prodR0 = velocitytep * deriv0_r_i; prodR1 = velocitytep * deriv1_r_i; prodR2 = velocitytep * deriv2_r_i;
magL0 = abs(prodL0); magL1 = abs(prodL1); magL2 = abs(prodL2);
magR0 = abs(prodR0); magR1 = abs(prodR1); magR2 = abs(prodR2);
flowL0 = calaFlowLValue(prodL0, prodR0, magL0, magR0);
flowL1 = calaFlowLValue(prodL1, prodR1, magL1, magR1);
flowL2 = calaFlowLValue(prodL2, prodR2, magL2, magR2);
flowR0 = calaFlowRValue(prodL0, prodR0, magL0, magR0);
flowR1 = calaFlowRValue(prodL1, prodR1, magL1, magR1);
flowR2 = calaFlowRValue(prodL2, prodR2, magL2, magR2);
//Now we know the upwind direction, add its contribution to \ | \grad \phi\ | .
FLOAT_TYPE magnitude =
pow(deriv0_l_i, 2) * flowL0 + pow(deriv0_r_i, 2) * flowR0 +
pow(deriv1_l_i, 2) * flowL1 + pow(deriv1_r_i, 2) * flowR1 +
pow(deriv2_l_i, 2) * flowL2 + pow(deriv2_r_i, 2) * flowR2;
// Finally, calculate speed* \ | \grad \phi\ |
magnitude = sqrtf(magnitude);
hamValue[i] = velocitytep * magnitude;
//magnitudeValue[begin_index + i] = magnitude;
FLOAT_TYPE effectiveVelocity0 = magL0 * flowL0 + magR0 * flowR0;
FLOAT_TYPE effectiveVelocity1 = magL1 * flowL1 + magR1 * flowR1;
FLOAT_TYPE effectiveVelocity2 = magL2 * flowL2 + magR2 * flowR2;
FLOAT_TYPE stepBoundInv0 = hji_grid->get_dxInv(0) * effectiveVelocity0;
FLOAT_TYPE stepBoundInv1 = hji_grid->get_dxInv(1) * effectiveVelocity1;
FLOAT_TYPE stepBoundInv2 = hji_grid->get_dxInv(2) * effectiveVelocity2;
if (magnitude > 0.0)
{
stepBoundMaxInv0 = stepBoundInv0 / magnitude > stepBoundMaxInv0 ? stepBoundInv0 / magnitude : stepBoundInv0;
stepBoundMaxInv1 = stepBoundInv1 / magnitude > stepBoundMaxInv1 ? stepBoundInv1 / magnitude : stepBoundInv1;
stepBoundMaxInv2 = stepBoundInv2 / magnitude > stepBoundMaxInv2 ? stepBoundInv2 / magnitude : stepBoundInv2;
}
}
step_bound_invs[0] = stepBoundMaxInv0;
step_bound_invs[1] = stepBoundMaxInv1;
step_bound_invs[2] = stepBoundMaxInv2;
return true;
}
void BurnNormalSchemeData::setVelocity(FLOAT_TYPE velocity)
{
this->velocity = velocity;
}
int BurnNormalSchemeData::calaFlowLValue(FLOAT_TYPE prodL, FLOAT_TYPE prodR, FLOAT_TYPE magL, FLOAT_TYPE magR) const
{
if ((prodL >= 0.0 && prodR >= 0.0) || (prodL >= 0.0 && prodR <= 0.0) && (magL >= magR))
{
return 1;
}
return 0;
}
int BurnNormalSchemeData::calaFlowRValue(FLOAT_TYPE prodL, FLOAT_TYPE prodR, FLOAT_TYPE magL, FLOAT_TYPE magR) const
{
if ((prodL <= 0.0 && prodR <= 0.0) || (prodL >= 0.0 && prodR <= 0.0) && (magL < magR))
{
return 1;
}
return 0;
}
- TermNormal.hpp
#ifndef __TermNormal_hpp__
#define __TermNormal_hpp__
//! Prefix to generate Visual C++ DLL
#ifdef _MAKE_VC_DLL
#define PREFIX_VC_DLL __declspec(dllexport)
//! Don't add prefix, except dll generating
#else
#define PREFIX_VC_DLL
#endif
#include <cstdint>
#include <vector>
#include <utility>
using namespace std::rel_ops;
#include <typedef.hpp>
#include <levelset/ExplicitIntegration/Terms/Term.hpp>
namespace levelset {
class SchemeData;
class TermNormal_impl;
class TermNormal:public Term
{
public:
/**
* motion in the normal direction in an HJ PDE with upwinding.
*
* @param [in] schemeData
* @param [in] type
* @return
*/
PREFIX_VC_DLL TermNormal(const SchemeData* schemeData, const beacls::UVecType type = beacls::UVecType_Vector);
/**
* @brief destructor
*
*/
~TermNormal();
/**
* .
*
* @param [out] ydot_ite
* @param [in] step_bound_invs
* @param [in] t
* @param [in] y
* @param [in] derivMins
* @param [in] derivMaxs
* @param [in] schemeData
* @param [in] loop_begin
* @param [in] loop_length
* @param [in] num_of_slices
* @param [in] enable_use_defined_dynamics_on_gpu
* @param [in] updataDerivMinMax
* @return true: Succeeded false :Failed Dissipation may be required global derivMins/derivMaxs.
Reduce partial derivMins/derivMaxs to global derivMins/derivMax, then execute again.
*/
bool execute(beacls::FloatVec::iterator ydot_ite,
beacls::FloatVec& step_bound_invs,
const FLOAT_TYPE t,
const beacls::FloatVec& y,
std::vector<beacls::FloatVec>& derivMins,
std::vector<beacls::FloatVec>& derivMaxs,
const SchemeData* schemeData,
const size_t loop_begin,
const size_t loop_length,
const size_t num_of_slices,
const bool enable_use_defined_dynamics_on_gpu,
const bool updataDerivMinMax
) const;
/**
* .
*
* @param schemeData
* @return true:Succeeded false:Failed
*/
bool synchronize(const SchemeData* schemeData) const;
bool operator==(const TermNormal& rhs) const;
bool operator==(const Term& rhs) const;
TermNormal* clone() const;
private:
TermNormal_impl* pimpl;
/**
* .
* @overload
* Disable operator =
* @param rhs
* @return
*/
TermNormal& operator=(const TermNormal& rhs);
/**
* .
* @overload
* Disable copy constructor
* @param rhs
*/
TermNormal(const TermNormal& rhs);
};
}
#endif /* __TermNormal_hpp__ */
###TermNormal_impl.hpp
#ifndef __TermNormal_impl_hpp__
#define __TermNormal_impl_hpp__
#include <cstdint>
#include <vector>
#include <Core/UVec.hpp>
#include <typedef.hpp>
namespace levelset {
class SchemeData;
class CacheTag;
class TermNormal_impl
{
private:
size_t first_dimension_loop_size;
size_t num_of_dimensions;
std::vector<beacls::UVec> deriv_l_uvecs;
std::vector<beacls::UVec> deriv_r_uvecs;
std::vector<beacls::UVec> deriv_c_uvecs;
std::vector<beacls::UVec> deriv_c_cpu_uvecs;
std::vector<beacls::UVec> x_uvecs;
//std::vector<beacls::UVec> deriv_max_uvecs;
//std::vector<beacls::UVec> deriv_min_uvecs;
beacls::UVec ydot_cuda_uvec;
/*maybe should delete*/
beacls::UVec ham_uvec;
beacls::UVec ham_cpu_uvec;
beacls::UVec diss_uvec;
/****************************************/
beacls::UVec flowR_uvec;
beacls::UVec flowL_uvec;
beacls::UVec magL_uvec;
beacls::UVec magR_uvec;
beacls::UVec magnitude_uvec;
beacls::UVecType type;
CacheTag* cacheTag;
public:
TermNormal_impl(
const SchemeData* schemeData,
const beacls::UVecType type);
~TermNormal_impl();
bool execute(
beacls::FloatVec::iterator ydot_ite,
beacls::FloatVec& step_bound_invs,
const FLOAT_TYPE t,
const beacls::FloatVec& y,
std::vector<beacls::FloatVec >& derivMins,
std::vector<beacls::FloatVec >& derivMaxs,
const SchemeData* schemeData,
const size_t loop_begin,
const size_t loop_length,
const size_t num_of_slices,
const bool enable_user_defined_dynamics_on_gpu,
const bool updateDerivMinMax
);
bool synchronize(const SchemeData* schemeData) const;
bool operator==(const TermNormal_impl& rhs) const {
if (this == &rhs) return true;
else if (type != rhs.type)return false;
else if (num_of_dimensions != rhs.num_of_dimensions)return false;
return true;
}
TermNormal_impl* clone() const {
return new TermNormal_impl(*this);
};
private:
/** @overload
Disable operator=
*/
TermNormal_impl& operator=(const TermNormal_impl& rhs);
/** @overload
Disable copy constructor
*/
TermNormal_impl(const TermNormal_impl& rhs);
};
#endif
}
- TermNormal.cpp
#include<levelset/ExplicitIntegration/Terms/TermNormal.hpp>
#include "TermNormal_impl.hpp"
#include <levelset/Grids/HJI_Grid.hpp>
//#include "TermNormal_cuda.hpp"//need to impl
#include <levelset/ExplicitIntegration/SchemeData.hpp>
#include <levelset/SpatialDerivative/SpatialDerivative.hpp>
#include <Core/CacheTag.hpp>
#include <algorithm>
#include <functional>
#include <typeinfo>
using namespace levelset;
using namespace std;
levelset::TermNormal_impl::TermNormal_impl(const SchemeData* schemeData, const beacls::UVecType type) :
first_dimension_loop_size(schemeData->get_grid()->get_N(0)),
num_of_dimensions(schemeData->get_grid()->get_num_of_dimensions()),
deriv_l_uvecs(num_of_dimensions),
deriv_r_uvecs(num_of_dimensions),
type(type),
cacheTag(new levelset::CacheTag())
{
}
levelset::TermNormal_impl::~TermNormal_impl()
{
if (cacheTag)
{
delete cacheTag;
}
}
levelset::TermNormal_impl::TermNormal_impl(const TermNormal_impl& rhs):
first_dimension_loop_size(rhs.first_dimension_loop_size),
num_of_dimensions(rhs.num_of_dimensions),
type(rhs.type),
cacheTag(new levelset::CacheTag())
{
deriv_l_uvecs.resize(rhs.deriv_l_uvecs.size());
deriv_r_uvecs.resize(rhs.deriv_r_uvecs.size());
x_uvecs.resize(rhs.x_uvecs.size());
}
bool levelset::TermNormal_impl::execute(
beacls::FloatVec::iterator ydot_ite,
beacls::FloatVec& step_bound_invs,
const FLOAT_TYPE t,
const beacls::FloatVec& y,
std::vector<beacls::FloatVec>& derivMins,
std::vector<beacls::FloatVec>& derivMaxs,
const SchemeData* schemeData,
const size_t loop_begin,
const size_t loop_length,
const size_t num_of_slices,
const bool enable_user_defined_dynamics_on_gpu,
const bool updateDerivMinMax)
{
const HJI_Grid* grid = schemeData->get_grid();
SpatialDerivative* spatialDerivative = schemeData->get_spatialDerivative();
beacls::UVecDepth depth = beacls::type_to_depth<FLOAT_TYPE>();
beacls::UVec y_uvec(y, beacls::UVecType_Vector, false);
size_t f_d_l_size = first_dimension_loop_size;//first dimension loop size
size_t grid_length = num_of_slices * loop_length * f_d_l_size;//
//size_t grid_length = grid->get_numel();
// printf("(num_of_slices, loop_length, f_d_l_size) = (%zu, %zu, %zu)\n",
// num_of_slices, loop_length, f_d_l_size);
// printf("grid_length = %zu\n", grid_length);
size_t slice_length = loop_length * f_d_l_size;
if (deriv_l_uvecs.size() != num_of_dimensions)
{
deriv_l_uvecs.resize(num_of_dimensions);
}
if (deriv_r_uvecs.size() != num_of_dimensions)
{
deriv_r_uvecs.resize(num_of_dimensions);
}
//init deriv value vector
for_each(deriv_l_uvecs.begin(), deriv_l_uvecs.end(), ([depth, grid_length, this](auto& rhs) {
if (rhs.type() != type)
{
rhs = beacls::UVec(depth, type, grid_length);
}
else if (rhs.size() != grid_length)
{
rhs.resize(grid_length);
}
}));
for_each(deriv_r_uvecs.begin(), deriv_r_uvecs.end(), ([depth, grid_length, this](auto& rhs) {
if (rhs.type() != type)
{
rhs = beacls::UVec(depth, type, grid_length);
}
else if (rhs.size() != grid_length)
{
rhs.resize(grid_length);
}
}));
//!< cala magnitude temp value
if (flowL_uvec.type() != type)
{
flowL_uvec = beacls::UVec(depth, type, grid_length);
}
else if (flowL_uvec.size() != grid_length)
{
flowL_uvec.resize(grid_length);
}
if (flowR_uvec.type() != type)
{
flowR_uvec = beacls::UVec(depth, type, grid_length);
}
else if (flowR_uvec.size() != grid_length)
{
flowR_uvec.resize(grid_length);
}
if (magL_uvec.type() != type)
{
magL_uvec = beacls::UVec(depth, type, grid_length);
}
else if (magL_uvec.size() != grid_length)
{
magL_uvec.resize(grid_length);
}
if (magR_uvec.type() != type)
{
magR_uvec = beacls::UVec(depth, type, grid_length);
}
else if (magR_uvec.size() != grid_length)
{
magR_uvec.resize(grid_length);
}
if (magnitude_uvec.type() != type)
{
magnitude_uvec = beacls::UVec(depth, type, grid_length);
}
else if (magnitude_uvec.size() != grid_length)
{
magnitude_uvec.resize(grid_length);
}
size_t src_index_term = loop_begin * f_d_l_size;
if (!cacheTag->check_tag(t, loop_begin, slice_length * num_of_slices))
{
//!< Copy xs to Cuda memory asynchronously in spatial derivative functions
x_uvecs.resize(num_of_dimensions);
//if use gpu init cuda vector
if (enable_user_defined_dynamics_on_gpu && (type == beacls::UVecType_Cuda))
{
for (size_t dimension = 0; dimension < num_of_dimensions; ++dimension) {
if (x_uvecs[dimension].type() != beacls::UVecType_Cuda) x_uvecs[dimension] = beacls::UVec(depth, beacls::UVecType_Cuda, grid_length);
else x_uvecs[dimension].resize(grid_length);
}
}
for (size_t index = 0; index < num_of_dimensions; ++index)
{
//!< To optimize asynchronous execution, calculate from heavy dimension (0, 2, 3 ... 1);
const size_t dimension = (index == 0) ? index : (index == num_of_dimensions - 1) ? 1 : index + 1;
beacls::UVec& deriv_l_uvec = deriv_l_uvecs[dimension];
beacls::UVec& deriv_r_uvec = deriv_r_uvecs[dimension];
spatialDerivative->execute(
deriv_l_uvec,
deriv_r_uvec,
grid,
y.data(),
dimension,
false,
loop_begin,
slice_length,
num_of_slices);
const beacls::FloatVec& xs = grid->get_xs(dimension);
beacls::copyHostPtrToUVecAsync(x_uvecs[dimension], xs.data() + src_index_term, grid_length);
}
for (size_t dimension = 0; dimension < num_of_dimensions; ++dimension)
{
beacls::UVec x_uvecs_dim = x_uvecs[dimension];
synchronizeUVec(x_uvecs_dim);
}
cacheTag->set_tag(t, loop_begin, slice_length * num_of_dimensions);
}
schemeData->normal3DHamFunc(
ham_uvec,
magnitude_uvec,
step_bound_invs,
deriv_l_uvecs,
deriv_r_uvecs,
src_index_term,
grid_length);
const beacls::FloatVec* ham_uvec_ptr = beacls::UVec_<FLOAT_TYPE>(ham_uvec).vec();
std::copy(ham_uvec_ptr->cbegin(), ham_uvec_ptr->cend(), ydot_ite);
return true;
}
bool levelset::TermNormal_impl::synchronize(const SchemeData* schemeData) const
{
beacls::synchronizeUVec(ydot_cuda_uvec);
return true;
}
levelset::TermNormal::TermNormal(const SchemeData* schemeData, const beacls::UVecType type)
{
pimpl = new TermNormal_impl(schemeData, type);
}
levelset::TermNormal::~TermNormal()
{
if (pimpl)
{
delete pimpl;
}
}
bool levelset::TermNormal::execute(
beacls::FloatVec::iterator ydot_ite,
beacls::FloatVec& step_bound_invs,
const FLOAT_TYPE t,
const beacls::FloatVec& y,
std::vector<beacls::FloatVec>& derivMins,
std::vector<beacls::FloatVec>& derivMaxs,
const SchemeData* schemeData,
const size_t loop_begin,
const size_t loop_length,
const size_t num_of_slices,
const bool enable_use_defined_dynamics_on_gpu,
const bool updataDerivMinMax) const
{
if (pimpl)
{
bool ret = pimpl->execute(ydot_ite, step_bound_invs, t, y, derivMins, derivMaxs, schemeData, loop_begin, loop_length, num_of_slices, enable_use_defined_dynamics_on_gpu, updataDerivMinMax);
return ret;
}
return false;
}
bool levelset::TermNormal::synchronize(const SchemeData* schemeData) const
{
if (pimpl)
{
return pimpl->synchronize(schemeData);
}
return false;
}
bool levelset::TermNormal::operator==(const TermNormal& rhs) const
{
if (this == &rhs) return true;
else if (!pimpl)
{
if (!rhs.pimpl)
{
return true;
}
else
{
return false;
}
}
else
{
if (!rhs.pimpl)
{
return false;
}
else
{
//object pointer address is equal
if (pimpl == rhs.pimpl)
{
return true;
}
else if (*pimpl == *rhs.pimpl)//object is equal
{
return true;
}
else
{
return false;
}
}
}
return false;
}
bool levelset::TermNormal::operator==(const Term& rhs) const
{
if (this == &rhs)
{
return true;
}
else if (typeid(*this) != typeid(rhs))//types are consistent
{
return false;
}
else
{
return operator==(dynamic_cast<const TermNormal&>(rhs));
}
}
TermNormal* levelset::TermNormal::clone() const
{
return new TermNormal(*this);
}
TermNormal& levelset::TermNormal::operator=(const TermNormal& rhs)
{
if (pimpl)
{
delete pimpl;
pimpl = rhs.pimpl->clone();
}
else
{
pimpl = rhs.pimpl->clone();
}
}
levelset::TermNormal::TermNormal(const TermNormal& rhs):pimpl(rhs.pimpl->clone())
{
}
- when i run this code and at 9.375s the data is start collapse as follow picture
YuDe95 commented
cala stepbound value cause this error. {use max{(stepBoundInv0+ stepBoundInv1+ stepBoundInv2) /mag}}