check grid values in delta eddington regression test

docker/Dockerfile.memcheck

@@ -26,6 +26,7 @@ RUN mkdir /build \
       && cd /build \
       && cmake \
             -DTUVX_ENABLE_MEMCHECK:BOOL=TRUE \
+            -DCMAKE_BUILD_TYPE=DEBUG \
         /tuv-x \
       && make -j 8
 

include/tuvx/util/array2d.hpp

@@ -34,6 +34,7 @@ namespace tuvx {
       typename std::vector<T>::const_iterator end() const { return data_.end(); }
 
       std::vector<T> &AsVector() { return data_; }
+      const std::vector<T> &AsVector() const { return data_; }
 
     private:
       size_t index(size_t i, size_t j) const { return i * dim2_ + j; }

include/tuvx/util/array3d.hpp

@@ -40,6 +40,7 @@ template <typename T = double> class Array3D {
   typename std::vector<T>::const_iterator end() const { return data_.end(); }
 
   std::vector<T> &AsVector() { return data_; }
+  const std::vector<T> &AsVector() const { return data_; }
 
 private:
   size_t index(size_t i, size_t j, size_t k) const {

test/regression/solvers/delta_eddington.F90

@@ -58,10 +58,6 @@ end subroutine free_output_c
   config_file_path = 'examples/tuv_5_4.json'
   call test_cpp_delta_eddington_solver_t(config_file_path)
 
-  ! Run the TS1/TSMLT test
-  config_file_path = 'examples/ts1_tsmlt.json'
-  call test_cpp_delta_eddington_solver_t(config_file_path)
-
   call musica_mpi_finalize( )
 
 contains
@@ -154,14 +150,14 @@ function calculate_cpp_radiation_fields( tuvx_core, solar_zenith_angles,      &
     solar_zenith_angles_c   = real( solar_zenith_angles(:), kind=c_double )   &
                                   * real(pi, kind=c_double) / 180.0_c_double ! degrees -> radians
     earth_sun_distances_c   = real( earth_sun_distances(:), kind=c_double )
-    allocate( altitude_mid_points_c( heights%ncells_,                         &
-                                     size( solar_zenith_angles ) ) )
-    allocate( altitude_edges_c(      heights%ncells_+1,                       &
-                                     size( solar_zenith_angles ) ) )
+    allocate( altitude_mid_points_c( size( solar_zenith_angles ),             &
+                                     heights%ncells_ ) )
+    allocate( altitude_edges_c(      size( solar_zenith_angles ),             &
+                                     heights%ncells_+1 ) )
     do i_column = 1, size( solar_zenith_angles )
-      altitude_mid_points_c(:,i_column) =                                     &
+      altitude_mid_points_c(i_column,:) =                                     &
           real( heights%mid_(:), kind=c_double ) * 1.0e3_c_double ! km -> m
-      altitude_edges_c(:,i_column)      =                                     &
+      altitude_edges_c(i_column,:)      =                                     &
           real( heights%edge_(:), kind=c_double ) * 1.0e3_c_double ! km -> m
     end do
     wavelength_mid_points_c = real( wavelengths%mid_(:), kind=c_double )      &
@@ -183,12 +179,12 @@ function calculate_cpp_radiation_fields( tuvx_core, solar_zenith_angles,      &
     do i_column = 1, size( solar_zenith_angles )
       n_lev = heights%ncells_
       n_wl = wavelengths%ncells_
-      allocate( radiation_fields(i_column)%edr_( n_lev, n_wl ) )
-      allocate( radiation_fields(i_column)%eup_( n_lev, n_wl ) )
-      allocate( radiation_fields(i_column)%edn_( n_lev, n_wl ) )
-      allocate( radiation_fields(i_column)%fdr_( n_lev, n_wl ) )
-      allocate( radiation_fields(i_column)%fup_( n_lev, n_wl ) )
-      allocate( radiation_fields(i_column)%fdn_( n_lev, n_wl ) )
+      allocate( radiation_fields(i_column)%edr_( n_lev+1, n_wl ) )
+      allocate( radiation_fields(i_column)%eup_( n_lev+1, n_wl ) )
+      allocate( radiation_fields(i_column)%edn_( n_lev+1, n_wl ) )
+      allocate( radiation_fields(i_column)%fdr_( n_lev+1, n_wl ) )
+      allocate( radiation_fields(i_column)%fup_( n_lev+1, n_wl ) )
+      allocate( radiation_fields(i_column)%fdn_( n_lev+1, n_wl ) )
       call copy_c_array_to_fortran( output%irrad_direct_,                      &
                                radiation_fields(i_column)%edr_,                &
                                i_column, heights%ncells_, wavelengths%ncells_ )

test/regression/solvers/delta_eddington.cpp

@@ -3,6 +3,12 @@
 //
 // Tests for tuvx::DeltaEddington.
 #include "delta_eddington.hpp"
+#include <iostream>
+
+#define ASSERT(x) if (!(x)) { \
+    std::cerr << "Assertion failed at " << __FILE__ << ":" << __LINE__ << " : "<< #x << std::endl; \
+    exit(EXIT_FAILURE); \
+}
 
 double* CopyVector(const std::vector<double>& vec)
 {
@@ -43,6 +49,49 @@ GridPolicy CreateFixedGrid(std::string units, const size_t sections, const doubl
   return grid;
 }
 
+// Checks certian input values to ensure the data is properly transferred
+// from the Fortran side to the C++ side.
+// If the testing conditions change in the future, this function will need to be updated.
+void CheckInputs(const std::vector<double>& solar_zenith_angles,
+                 const std::vector<double>& earth_sun_distances,
+                 const std::map<std::string, tuvx::Grid<tuvx::Array2D<double>>>& grids,
+                 const std::map<std::string, tuvx::Profile<tuvx::Array2D<double>>>& profiles,
+                 const tuvx::RadiatorState<tuvx::Array3D<double>>& accumulated_radiator_states)
+{
+  ASSERT(solar_zenith_angles.size() == 2);
+  ASSERT(earth_sun_distances.size() == 2);
+  ASSERT(grids.size() == 2);
+  ASSERT(grids.at("altitude [m]").NumberOfColumns() == 2);
+  ASSERT(grids.at("altitude [m]").NumberOfSections() == 120);
+  ASSERT(!grids.at("altitude [m]").IsConstant());
+  // heights have been converted to meters from km
+  ASSERT(grids.at("altitude [m]").edges_(0,0) == 0.0); 
+  ASSERT(grids.at("altitude [m]").edges_(42,0) == 42000.0);
+  ASSERT(grids.at("altitude [m]").edges_(120,0) == 120000.0); 
+  ASSERT(grids.at("altitude [m]").edges_(0,1) == 0.0); 
+  ASSERT(grids.at("altitude [m]").edges_(42,1) == 42000.0);
+  ASSERT(grids.at("altitude [m]").edges_(120,1) == 120000.0); 
+  ASSERT(grids.at("altitude [m]").mid_points_(0,0) == 500.0);
+  ASSERT(grids.at("altitude [m]").mid_points_(41,0) == 41500.0);
+  ASSERT(grids.at("altitude [m]").mid_points_(119,0) == 119500.0);
+  ASSERT(grids.at("altitude [m]").mid_points_(0,1) == 500.0);
+  ASSERT(grids.at("altitude [m]").mid_points_(41,1) == 41500.0);
+  ASSERT(grids.at("altitude [m]").mid_points_(119,1) == 119500.0);
+  ASSERT(grids.at("wavelength [m]").NumberOfColumns() == 1);
+  ASSERT(grids.at("wavelength [m]").NumberOfSections() == 156);
+  ASSERT(grids.at("wavelength [m]").IsConstant());
+  // wavelengths have been converted to meters from nm
+  ASSERT(grids.at("wavelength [m]").edges_(0,0) == 120.0*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").edges_(77,0) == 311.5*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").edges_(156,0) == 735.0*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").mid_points_(0,0) > 120.69*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").mid_points_(0,0) < 120.71*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").mid_points_(77,0) > 311.9*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").mid_points_(77,0) < 312.1*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").mid_points_(155,0) > 729.9*1.0e-9);
+  ASSERT(grids.at("wavelength [m]").mid_points_(155,0) < 730.1*1.0e-9);
+}
+
 SolverOutput RunDeltaEddingtonSolver(const SolverInput input)
 {
   using GridPolicy = tuvx::Grid<tuvx::Array2D<double>>;
@@ -59,6 +108,7 @@ SolverOutput RunDeltaEddingtonSolver(const SolverInput input)
                                                                                               grids["altitude [m]"],
                                                                                               grids["wavelength [m]"]);
   tuvx::DeltaEddington solver;
+  CheckInputs(solar_zenith_angles, earth_sun_distances, grids, profiles, accumulated_radiator_states);
   solver.Solve(solar_zenith_angles, grids, profiles, accumulated_radiator_states, radiation_field);
   SolverOutput output;
   output.n_wavelengths_ = input.n_wavelengths_;