Fix netcdf output bug for timestep and timeblock. closes #257 (#333)

jcurtis2 · web-flow · commit 6b59486ccd8e · 2024-01-31T20:08:22.000Z
diff --git a/examples/mie_optical.ipynb b/examples/mie_optical.ipynb
diff --git a/examples/particle_simulation.ipynb b/examples/particle_simulation.ipynb
diff --git a/src/run_part.F90 b/src/run_part.F90
@@ -86,7 +86,10 @@ subroutine f_run_part_timestep( &
     camp_core_ptr_c, &
     photolysis_ptr_c, &
     i_time, &
-    t_start &
+    t_start, &
+    last_output_time, &
+    last_progress_time, &
+    i_output &
   ) bind(C)
 
     type(c_ptr), intent(in) :: scenario_ptr_c
@@ -117,8 +120,11 @@ subroutine f_run_part_timestep( &
     integer(c_int), intent(in) :: i_time
     real(c_double), intent(in) :: t_start
 
-    real(c_double) :: last_output_time, last_progress_time
-    integer(c_int) :: i_output, progress_n_samp, progress_n_coag, &
+    real(c_double), intent(inout) :: last_output_time
+    real(c_double), intent(inout) :: last_progress_time
+    integer(c_int), intent(inout) :: i_output
+
+    integer(c_int) :: progress_n_samp, progress_n_coag, &
         progress_n_emit, progress_n_dil_in, progress_n_dil_out, &
         progress_n_nuc
 
@@ -142,6 +148,13 @@ subroutine f_run_part_timestep( &
     if (env_state_ptr_f%elapsed_time < run_part_opt_ptr_f%del_t) then
        call mosaic_init(env_state_ptr_f, aero_data_ptr_f, run_part_opt_ptr_f%del_t, &
             run_part_opt_ptr_f%do_optical)
+       if (run_part_opt_ptr_f%t_output > 0) then
+          call output_state(run_part_opt_ptr_f%output_prefix, &
+               run_part_opt_ptr_f%output_type, aero_data_ptr_f, aero_state_ptr_f, gas_data_ptr_f, &
+               gas_state_ptr_f, env_state_ptr_f, 1, .0d0, run_part_opt_ptr_f%del_t, &
+               run_part_opt_ptr_f%i_repeat, run_part_opt_ptr_f%record_removals, &
+               run_part_opt_ptr_f%do_optical, run_part_opt_ptr_f%uuid)
+       end if
     end if
     call run_part_timestep(scenario_ptr_f, env_state_ptr_f, aero_data_ptr_f, aero_state_ptr_f, &
        gas_data_ptr_f, gas_state_ptr_f, run_part_opt_ptr_f, camp_core_ptr_f, photolysis_ptr_f, &
@@ -164,7 +177,10 @@ subroutine f_run_part_timeblock( &
     photolysis_ptr_c, &
     i_time, &
     i_next, &
-    t_start &
+    t_start, &
+    last_output_time, &
+    last_progress_time, &
+    i_output &
   ) bind(C)
 
     type(c_ptr), intent(in) :: scenario_ptr_c
@@ -197,9 +213,11 @@ subroutine f_run_part_timeblock( &
     integer(c_int), intent(in) :: i_time
     integer(c_int), intent(in) :: i_next
     real(c_double), intent(in) :: t_start
+    real(c_double), intent(inout) :: last_output_time
+    real(c_double), intent(inout) :: last_progress_time
+    integer(c_int), intent(inout) :: i_output
 
-    real(c_double) :: last_output_time, last_progress_time
-    integer(c_int) :: i_output, progress_n_samp, progress_n_coag, &
+    integer(c_int) :: progress_n_samp, progress_n_coag, &
         progress_n_emit, progress_n_dil_in, progress_n_dil_out, &
         progress_n_nuc
 
@@ -223,6 +241,13 @@ subroutine f_run_part_timeblock( &
     if (env_state_ptr_f%elapsed_time < run_part_opt_ptr_f%del_t) then
        call mosaic_init(env_state_ptr_f, aero_data_ptr_f, run_part_opt_ptr_f%del_t, &
             run_part_opt_ptr_f%do_optical)
+       if (run_part_opt_ptr_f%t_output > 0) then
+          call output_state(run_part_opt_ptr_f%output_prefix, &
+               run_part_opt_ptr_f%output_type, aero_data_ptr_f, aero_state_ptr_f, gas_data_ptr_f, &
+               gas_state_ptr_f, env_state_ptr_f, 1, .0d0, run_part_opt_ptr_f%del_t, &
+               run_part_opt_ptr_f%i_repeat, run_part_opt_ptr_f%record_removals, &
+               run_part_opt_ptr_f%do_optical, run_part_opt_ptr_f%uuid)
+       end if
     end if
 
     call run_part_timeblock(scenario_ptr_f, env_state_ptr_f, aero_data_ptr_f, aero_state_ptr_f, &
diff --git a/src/run_part.cpp b/src/run_part.cpp
@@ -5,6 +5,7 @@
 ##################################################################################################*/
 
 #include "run_part.hpp"
+#include "pybind11/stl.h"
 
 void run_part(
     const Scenario &scenario,
@@ -30,7 +31,7 @@ void run_part(
     );
 }
 
-void run_part_timestep(
+std::tuple<double, double, int> run_part_timestep(
     const Scenario &scenario,
     EnvState &env_state,
     const AeroData &aero_data,
@@ -41,7 +42,10 @@ void run_part_timestep(
     const CampCore &camp_core,
     const Photolysis &photolysis,
     const int &i_time,
-    const double &t_start
+    const double &t_start,
+    double &last_output_time,
+    double &last_progress_time,
+    int &i_output
 ) {
     f_run_part_timestep(
         scenario.ptr.f_arg(),
@@ -54,11 +58,16 @@ void run_part_timestep(
         camp_core.ptr.f_arg(),
         photolysis.ptr.f_arg(),
         &i_time,
-        &t_start
+        &t_start,
+        &last_output_time,
+        &last_progress_time,
+        &i_output
     );
+
+    return std::make_tuple(last_output_time, last_progress_time, i_output);
 }
 
-void run_part_timeblock(
+std::tuple<double, double, int> run_part_timeblock(
     const Scenario &scenario,
     EnvState &env_state,
     const AeroData &aero_data,
@@ -70,7 +79,10 @@ void run_part_timeblock(
     const Photolysis &photolysis,
     const int &i_time,
     const int &i_next,
-    const double &t_start
+    const double &t_start,
+    double &last_output_time,
+    double &last_progress_time,
+    int &i_output
 ) {
     f_run_part_timeblock(
         scenario.ptr.f_arg(),
@@ -84,6 +96,11 @@ void run_part_timeblock(
         photolysis.ptr.f_arg(),
         &i_time,
         &i_next,
-        &t_start
+        &t_start,
+        &last_output_time,
+        &last_progress_time,
+        &i_output
     );
+
+    return std::make_tuple(last_output_time, last_progress_time, i_output);
 }
diff --git a/src/run_part.hpp b/src/run_part.hpp
@@ -38,7 +38,10 @@ extern "C" void f_run_part_timestep(
     const void*,
     const void*,
     const int *,
-    const double *
+    const double*,
+    double*,
+    double*,
+    int*
 ) noexcept;
 
 extern "C" void f_run_part_timeblock(
@@ -53,7 +56,10 @@ extern "C" void f_run_part_timeblock(
     const void*,
     const int *,
     const int *,
-    const double *
+    const double *,
+    double*,
+    double*,
+    int*
 ) noexcept;
 
 void run_part(
@@ -68,7 +74,7 @@ void run_part(
     const Photolysis &photolysis
 );
 
-void run_part_timestep(
+std::tuple<double, double, int> run_part_timestep(
     const Scenario &scenario,
     EnvState &env_state,
     const AeroData &aero_data,
@@ -79,10 +85,13 @@ void run_part_timestep(
     const CampCore &camp_core,
     const Photolysis &photolysis,
     const int &i_time,
-    const double &t_start
+    const double &t_start,
+    double &last_output_time,
+    double &last_progress_time,
+    int &i_output
 );
 
-void run_part_timeblock(
+std::tuple<double, double, int> run_part_timeblock(
     const Scenario &scenario,
     EnvState &env_state,
     const AeroData &aero_data,
@@ -94,5 +103,8 @@ void run_part_timeblock(
     const Photolysis &photolysis,
     const int &i_time,
     const int &i_next,
-    const double &t_start
+    const double &t_start,
+    double &last_output_time,
+    double &last_progress_time,
+    int &i_output
 );
diff --git a/tests/test_run_part.py b/tests/test_run_part.py
@@ -12,26 +12,30 @@
 from .test_aero_state import AERO_STATE_CTOR_ARG_MINIMAL
 from .test_env_state import ENV_STATE_CTOR_ARG_MINIMAL
 from .test_gas_data import GAS_DATA_CTOR_ARG_MINIMAL
-from .test_run_part_opt import RUN_PART_OPT_CTOR_ARG_MINIMAL
+from .test_run_part_opt import RUN_PART_OPT_CTOR_ARG_SIMULATION
 from .test_scenario import SCENARIO_CTOR_ARG_MINIMAL
 
 
 @pytest.fixture
-def common_args():
+def common_args(tmp_path):
     aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
     gas_data = ppmc.GasData(GAS_DATA_CTOR_ARG_MINIMAL)
     gas_state = ppmc.GasState(gas_data)
     scenario = ppmc.Scenario(gas_data, aero_data, SCENARIO_CTOR_ARG_MINIMAL)
     env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL)
     scenario.init_env_state(env_state, 0.0)
+    filename = tmp_path / "test"
+    run_part_opt = ppmc.RunPartOpt(
+        {**RUN_PART_OPT_CTOR_ARG_SIMULATION, "output_prefix": str(filename)}
+    )
     return (
         scenario,
         env_state,
         aero_data,
         ppmc.AeroState(aero_data, *AERO_STATE_CTOR_ARG_MINIMAL),
         gas_data,
         gas_state,
-        ppmc.RunPartOpt(RUN_PART_OPT_CTOR_ARG_MINIMAL),
+        run_part_opt,
         ppmc.CampCore(),
         ppmc.Photolysis(),
     )
@@ -42,10 +46,29 @@ class TestRunPart:
     def test_run_part(common_args):  # pylint: disable=redefined-outer-name
         ppmc.run_part(*common_args)
 
+        assert common_args[1].elapsed_time == RUN_PART_OPT_CTOR_ARG_SIMULATION["t_max"]
+
     @staticmethod
     def test_run_part_timestep(common_args):  # pylint: disable=redefined-outer-name
-        ppmc.run_part_timestep(*common_args, 0, 0)
+        (last_output_time, last_progress_time, i_output) = ppmc.run_part_timestep(
+            *common_args, 1, 0, 0, 0, 1
+        )
+
+        assert common_args[1].elapsed_time == RUN_PART_OPT_CTOR_ARG_SIMULATION["del_t"]
+        assert last_output_time == 0.0
+        assert last_progress_time == 0.0
+        assert i_output == 1
 
     @staticmethod
     def test_run_part_timeblock(common_args):  # pylint: disable=redefined-outer-name
-        ppmc.run_part_timeblock(*common_args, 0, 0, 0)
+
+        num_times = int(
+            RUN_PART_OPT_CTOR_ARG_SIMULATION["t_output"]
+            / RUN_PART_OPT_CTOR_ARG_SIMULATION["del_t"]
+        )
+        (last_output_time, last_progress_time, i_output) = ppmc.run_part_timeblock(
+            *common_args, 1, num_times, 0, 0, 0, 1
+        )
+        assert last_output_time == RUN_PART_OPT_CTOR_ARG_SIMULATION["t_output"]
+        assert last_progress_time == 0.0
+        assert i_output == 2
