Merge pull request #471 from arkedge/feature/add_thermal2

 Solve fixmes and todos in thermal simulation (Re-re-open, change base repository)

src/dynamics/dynamics.cpp

@@ -29,7 +29,8 @@ void Dynamics::Initialize(const SimulationConfiguration* simulation_configuratio
                      local_celestial_information.GetGlobalInformation().GetCenterBodyGravityConstant_m3_s2(), "ORBIT", relative_information);
   attitude_ = InitAttitude(simulation_configuration->spacecraft_file_list_[spacecraft_id], orbit_, &local_celestial_information,
                            simulation_time->GetAttitudeRkStepTime_s(), structure->GetKinematicsParameters().GetInertiaTensor_b_kgm2(), spacecraft_id);
-  temperature_ = InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s());
+  temperature_ = InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s(),
+                                 &(local_environment_->GetSolarRadiationPressure()));
 
   // To get initial value
   orbit_->UpdateByAttitude(attitude_->GetQuaternion_i2b());

src/dynamics/thermal/heater.cpp

@@ -1,10 +1,12 @@
 #include "heater.hpp"
 
+#include <cassert>
 #include <cmath>
 
 using namespace std;
 
 Heater::Heater(const unsigned int heater_id, const double power_rating_W) : heater_id_(heater_id), power_rating_W_(power_rating_W) {
+  AssertHeaterParams();
   heater_status_ = HeaterStatus::kOff;
   power_output_W_ = 0;
 }
@@ -26,3 +28,8 @@ void Heater::PrintParam(void) {
   cout << "  status            : " << static_cast<int>(heater_status_) << endl;
   cout << "  power output      : " << power_output_W_ << endl;
 }
+
+void Heater::AssertHeaterParams(void) {
+  assert(heater_id_ >= 1);       // Heater ID must be larger than 1
+  assert(power_rating_W_ >= 0);  // Power rating must be larger than 0
+}

src/dynamics/thermal/heater.hpp

@@ -83,6 +83,12 @@ class Heater {
    * @brief Print power_rating_W, heater_status_, power_output_W_
    */
   void PrintParam(void);
+
+  /**
+   * @fn AssertHeaterParams
+   * @brief Check Heater Parameters
+   */
+  void AssertHeaterParams(void);
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_HEATER_HPP_

src/dynamics/thermal/heater_controller.cpp

@@ -1,11 +1,14 @@
 #include "heater_controller.hpp"
 
+#include <cassert>
 #include <cmath>
 
 using namespace std;
 
 HeaterController::HeaterController(const double lower_threshold_degC, const double upper_threshold_degC)
-    : lower_threshold_degC_(lower_threshold_degC), upper_threshold_degC_(upper_threshold_degC) {}
+    : lower_threshold_degC_(lower_threshold_degC), upper_threshold_degC_(upper_threshold_degC) {
+  AssertHeaterControllerParams();
+}
 
 HeaterController::~HeaterController() {}
 
@@ -21,3 +24,5 @@ void HeaterController::ControlHeater(Heater* p_heater, double temperature_degC)
     }
   }
 }
+
+void HeaterController::AssertHeaterControllerParams() { assert(upper_threshold_degC_ > lower_threshold_degC_); }

src/dynamics/thermal/heater_controller.hpp

@@ -10,6 +10,7 @@
 #include <library/logger/logger.hpp>
 #include <string>
 #include <vector>
+
 #include "heater.hpp"
 
 class HeaterController {
@@ -72,6 +73,12 @@ class HeaterController {
    * @param[in] temperature_degC
    */
   void ControlHeater(Heater* p_heater, double temperature_degC);
+
+  /**
+   * @fn AssertHeaterControllerParams
+   * @brief Check Parameters of HeaterController
+   */
+  void AssertHeaterControllerParams(void);
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_HEATER_CONTROLLER_HPP_

src/dynamics/thermal/heatload.cpp

@@ -1,5 +1,6 @@
 #include "heatload.hpp"
 
+#include <cassert>
 #include <cmath>
 #include <environment/global/physical_constants.hpp>
 
@@ -15,8 +16,8 @@ Heatload::Heatload(int node_id, std::vector<double> time_table_s, std::vector<do
   heater_heatload_W_ = 0.0;
   total_heatload_W_ = 0.0;
 
-  // [FIXME] Include assertion for size of time_table_s_ and internal_heatload_table_W_ to be larger than 2
-  // [FIXME] Include assertion for size of time_table_s_ and internal_heatload_table_W_ to be the same
+  AssertHeatloadParams();
+
   time_table_period_s_ = time_table_s_[time_table_s_.size() - 1];
   residual_elapsed_time_s_ = fmod(elapsed_time_s_, time_table_period_s_);
 }
@@ -48,3 +49,11 @@ void Heatload::SetElapsedTime_s(double elapsed_time_s) {
     }
   }
 }
+
+void Heatload::AssertHeatloadParams() {
+  // size of time_table_s_ and internal_heatload_table_W_ must be larger than 1
+  assert(time_table_s_.size() >= 1);
+  assert(internal_heatload_table_W_.size() >= 1);
+  // size of time_table_s_ and internal_heatload_table_W_ must be the same
+  assert(time_table_s_.size() == internal_heatload_table_W_.size());
+}

src/dynamics/thermal/heatload.hpp

@@ -99,6 +99,12 @@ class Heatload {
    * @param[in] heater_heatload_W
    */
   inline void SetHeaterHeatload_W(double heater_heatload_W) { heater_heatload_W_ = heater_heatload_W; }
+
+  /**
+   * @fn AssertHeatloadParams
+   * @brief Check if Heatload Parameters are Correct
+   */
+  void AssertHeatloadParams();
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_HEATLOAD_HPP_

src/dynamics/thermal/initialize_heater.cpp

@@ -5,12 +5,26 @@
 
 #include "initialize_heater.hpp"
 
-// [FIXME] Write discription of csv file for heater initialization
+#include <cassert>
+
+/* Import heater properties by reading CSV File (heaters.csv)
+
+[File Formats of heater.csv]
+column 1: Heater_id(int, Use values larger than or equal to 1)
+column 2: Power Rating (double, [W])
+column 3: Lower threshold of control (double, [degC])
+column 4: Upper threshold of control (double, [degC])
+
+First row is for Header, data begins from the second row
+*/
 
 Heater InitHeater(const std::vector<std::string>& heater_str) {
   using std::stod;
   using std::stoi;
 
+  size_t heater_str_size_defined = 4;                    // Correct size of heater_str
+  assert(heater_str.size() == heater_str_size_defined);  // Check if size of heater_str is correct
+
   int heater_id = 0;
   double power_rating_W = 0;  // [W]
 
@@ -28,6 +42,9 @@ Heater InitHeater(const std::vector<std::string>& heater_str) {
 HeaterController InitHeaterController(const std::vector<std::string>& heater_str) {
   using std::stod;
 
+  size_t heater_str_size_defined = 4;                    // Correct size of heater_str
+  assert(heater_str.size() == heater_str_size_defined);  // Check if size of heater_str is correct
+
   // Index to read from heater_str for each parameter
   int index_lower_threshold = 2;
   int index_upper_threshold = 3;

src/dynamics/thermal/initialize_heatload.cpp

@@ -5,17 +5,40 @@
 
 #include "initialize_heatload.hpp"
 
-// [FIXME] Write discription of csv file for heatload initialization
+#include <cassert>
+
+/* Import heatload properties by reading CSV File (heatload.csv)
+
+[File Formats of heatload.csv]
+First Row: Array of Simulation Times [s]
+After Second Row:
+  First Column: Node ID
+  After Second Column: Array of Heatload Values [W]
+
+Heatload object calculates heatload value of certain time by interpolating values of the heatload array
+
+Data Example:
+NodeID/Time, 0, 99, 100, 300, 500
+0,          10, 10,  20,  10, 20
+
+In this example, heatload at time 200s is calculated as 15W from interpolation
+
+The number of nodes and heatloads must be the same
+
+*/
 
 Heatload InitHeatload(const std::vector<std::string>& time_str, const std::vector<std::string>& internal_heatload_str) {
   using std::stod;
   using std::stoi;
 
-  // [FIXME] Include assertion for size of time_str and internal_heatload_str to be larger than 2
+  assert(time_str.size() >= 2);                             // Number of columns must be larger than 2
+  assert(internal_heatload_str.size() >= 2);                // Number of columns must be larger than 2
+  assert(time_str.size() == internal_heatload_str.size());  // Number of columns must be same
+
   std::vector<double> time_table_s(time_str.size() - 1);                            // exclude index
   std::vector<double> internal_heatload_table_W(internal_heatload_str.size() - 1);  // exclude index
 
-  int node_id = stoi(internal_heatload_str[0]);                                     // First data of internal_heatload_str is node id
+  int node_id = stoi(internal_heatload_str[0]);  // First data of internal_heatload_str is node id
 
   // read table
   for (unsigned int i = 0; i < time_str.size() - 1; ++i) {

src/dynamics/thermal/initialize_node.cpp

@@ -5,15 +5,16 @@
 
 #include "initialize_node.hpp"
 
+#include <cassert>
 #include <iostream>
 #include <library/initialize/initialize_file_access.hpp>
 #include <string>
 #include <typeinfo>
 #include <vector>
 
-/* Import node properties by reading CSV File (Node.csv)
+/* Import node properties by reading CSV File (node.csv)
 
-[File Formats of Node.csv]
+[File Formats of node.csv]
 column 1: Node_id(int)
 column 2: Node_label(string)
 column 3: Node_type(int, 0: diffusive, 1: boundary), Arithmetic node to be implemented as future work
@@ -38,6 +39,9 @@ Node InitNode(const std::vector<std::string>& node_str) {
   using std::stod;
   using std::stoi;
 
+  size_t node_str_size_defined = 11;                 // Correct size of node_str
+  assert(node_str.size() == node_str_size_defined);  // Check if size of node_str is correct
+
   int node_id = 0;                  // node number
   std::string node_label = "temp";  // node name
   int node_type_int = 0;            // node type
@@ -69,6 +73,15 @@ Node InitNode(const std::vector<std::string>& node_str) {
   for (int i = 0; i < 3; i++) {
     normal_v_b[i] = stod(node_str[index_normal_v_b_head + i]);
   }
+
+  // Normalize Norm Vector (Except for Boundary and Arithmetic Nodes)
+  if (node_type_int == 0) {
+    double norm = normal_v_b.CalcNorm();
+    for (int i = 0; i < 3; i++) {
+      normal_v_b[i] = normal_v_b[i] / norm;
+    }
+  }
+
   temperature_K = stod(node_str[index_temperature]);
 
   NodeType node_type = NodeType::kDiffusive;

src/dynamics/thermal/initialize_temperature.cpp

@@ -5,6 +5,7 @@
 
 #include "initialize_temperature.hpp"
 
+#include <cassert>
 #include <environment/global/simulation_time.hpp>
 #include <library/initialize/initialize_file_access.hpp>
 #include <string>
@@ -57,7 +58,7 @@ First row is time data
 using std::string;
 using std::vector;
 
-Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s) {
+Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s, const SolarRadiationPressureEnvironment* srp_environment) {
   auto mainIni = IniAccess(file_name);
 
   vector<Node> node_list;
@@ -70,6 +71,7 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
   vector<vector<string>> heater_str_list;    // string vector of heater property data
   vector<vector<string>> heatload_str_list;  // string vector of heatload property data
   size_t node_num = 1;
+  size_t heatload_num = 1;
   size_t heater_num = 1;
 
   bool is_calc_enabled = mainIni.ReadEnable("THERMAL", "calculation");
@@ -100,6 +102,7 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
   /*since we don't know the number of node_list yet, set node_num=100 temporary.
     Recall that Nodes_num are given to this function only to reseve memory*/
 
+  heatload_num = heatload_str_list.size() - 1;
   auto times_itr = heatload_str_list.begin();  // First Row is Time Data
   for (auto itr = heatload_str_list.begin() + 1; itr != heatload_str_list.end(); ++itr) {
     heatload_list.push_back(InitHeatload(*times_itr, *itr));
@@ -118,6 +121,8 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
     node_list.push_back(InitNode(*itr));
   }
 
+  assert(node_num == heatload_num);  // Number of nodes and heatload lists must be the same
+
   // Read Heater Properties from CSV File
   string filepath_heater = file_path + "heaters.csv";
   IniAccess conf_heater(filepath_heater);
@@ -141,8 +146,13 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
   conf_cij.ReadCsvDoubleWithHeader(conductance_matrix, node_num, 1, 1);
   conf_rij.ReadCsvDoubleWithHeader(radiation_matrix, node_num, 1, 1);
 
+  assert(conductance_matrix.size() == node_num);                      // Dimension must be same as node_num
+  assert(radiation_matrix.size() == node_num);                        // Dimension must be same as node_num
+  assert(conductance_matrix.size() == conductance_matrix[0].size());  // Must be square matrix
+  assert(radiation_matrix.size() == radiation_matrix[0].size());      // Must be square matrix
+
   Temperature* temperature;
   temperature = new Temperature(conductance_matrix, radiation_matrix, node_list, heatload_list, heater_list, heater_controller_list, node_num,
-                                rk_prop_step_s, is_calc_enabled, solar_calc_setting, debug);
+                                rk_prop_step_s, srp_environment, is_calc_enabled, solar_calc_setting, debug);
   return temperature;
 }

src/dynamics/thermal/initialize_temperature.hpp

@@ -15,6 +15,6 @@
  * @param[in] rk_prop_step_s: time step interval for temperature propagation integration
  * @return Temperature*
  */
-Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s);
+Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s, const SolarRadiationPressureEnvironment* srp_environment);
 
 #endif  // S2E_DYNAMICS_THERMAL_INITIALIZE_TEMPERATURE_HPP_

src/dynamics/thermal/node.cpp

@@ -5,6 +5,7 @@
 
 #include "node.hpp"
 
+#include <cassert>
 #include <cmath>
 #include <environment/global/physical_constants.hpp>
 
@@ -20,25 +21,20 @@ Node::Node(const int node_id, const string node_name, const NodeType node_type,
       capacity_J_K_(capacity_J_K),
       alpha_(alpha),
       area_m2_(area_m2),
-      normal_vector_b_(normal_vector_b),
-      node_type_(node_type) {
+      node_type_(node_type),
+      normal_vector_b_(normal_vector_b) {
+  AssertNodeParams();
   solar_radiation_W_ = 0;
 }
 
 Node::~Node() {}
 
-double Node::CalcSolarRadiation_W(libra::Vector<3> sun_direction_b) {
-  // FIXME: constants
-  double R = 6.96E+8;                              // Distance from sun
-  double T = 5778;                                 // sun surface temperature [K]
-  double sigma = 5.67E-8;                          // stephan-boltzman constant
-  double a = sun_direction_b.CalcNorm();           // distance from sun
-  double S = pow((R / a), 2) * sigma * pow(T, 4);  // Solar Constant at s/c position S[W/m^2]
-  double cos_theta = InnerProduct(sun_direction_b, normal_vector_b_) / a / normal_vector_b_.CalcNorm();
+double Node::CalcSolarRadiation_W(libra::Vector<3> sun_direction_b, double solar_flux_W_m2) {
+  double cos_theta = InnerProduct(sun_direction_b, normal_vector_b_);
 
   // calculate sun_power
   if (cos_theta > 0)
-    solar_radiation_W_ = S * area_m2_ * alpha_ * cos_theta;
+    solar_radiation_W_ = solar_flux_W_m2 * area_m2_ * alpha_ * cos_theta;
   else
     solar_radiation_W_ = 0;
   return solar_radiation_W_;
@@ -66,3 +62,12 @@ void Node::PrintParam(void) {
   }
   cout << endl;
 }
+
+void Node::AssertNodeParams(void) {
+  assert(node_id_ >= 0);                                      // Node ID should be larger than 0
+  assert(heater_id_ >= 0);                                    // Heater ID should be larger than 0
+  assert(temperature_K_ >= environment::absolute_zero_degC);  // Temperature must be larger than zero kelvin
+  assert(capacity_J_K_ >= 0);                                 // Capacity must be larger than 0, use 0 when node is boundary or arithmetic
+  assert(0 <= alpha_ && alpha_ <= 1);                         // alpha must be between 0 and 1
+  assert(area_m2_ >= 0);                                      // Area must be larger than 0
+}

src/dynamics/thermal/node.hpp

@@ -70,7 +70,7 @@ class Node {
    * @param sun_direction_b: Sun direction in body frame
    * @return double: Solar Radiation [W]
    */
-  double CalcSolarRadiation_W(libra::Vector<3> sun_direction_b);  // 太陽入射熱を計算
+  double CalcSolarRadiation_W(libra::Vector<3> sun_direction_b, double solar_flux_W_m2);
 
   // Getter
   /**
@@ -137,6 +137,12 @@ class Node {
    * @brief Print parameters of node in debug output
    */
   void PrintParam(void);
+
+  /**
+   * @fn AssertNodeParams
+   * @brief Check if Node Parameters are Correct
+   */
+  void AssertNodeParams(void);
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_NODE_HPP_