Framework and Executable Changes (Code Cleanup)

examples/CMakeLists.txt

@@ -1,6 +1,4 @@
-add_executable(PFHub1a_Periodic PFHub1a_Periodic.cpp)
-target_link_libraries(PFHub1a_Periodic LINK_PUBLIC CabanaPF)
-add_executable(PFHub1a_Benchmark PFHub1a_Benchmark.cpp)
-target_link_libraries(PFHub1a_Benchmark LINK_PUBLIC CabanaPF)
+add_executable(PFHub1a PFHub1a.cpp)
+target_link_libraries(PFHub1a LINK_PUBLIC CabanaPF)
 add_executable(PerformanceRuns GridTimer.cpp)
 target_link_libraries(PerformanceRuns LINK_PUBLIC CabanaPF)

examples/GridTimer.cpp

@@ -12,12 +12,17 @@ int main(int argc, char* argv[]) {
 
         // read in grid points to run off the command line:
         std::vector<int> runs;
+        double dt = 0, end_time = 0;
         try {
-            for (int i = 1; i < argc; i++)
+            if (argc < 4)
+                // Need at least: Executable name, timesteps, end time, one grid points
+                throw std::invalid_argument("");
+            dt = std::stoi(argv[2]);
+            end_time = std::stoi(argv[3]);
+            for (int i = 3; i < argc; i++)
                 runs.push_back(std::stoi(argv[i]));
-        } catch (std::logic_error const&) {
-            std::cout << "Usage: ./GridTimer [grid_points] [grid_points] [...]" << std::endl;
-            ;
+        } catch (std::invalid_argument const&) {
+            std::cout << "Usage: ./PerformanceRuns <dt> <end time> <grid points> [grid points ...]" << std::endl;
             return 1;
         }
 
@@ -27,8 +32,8 @@ int main(int argc, char* argv[]) {
             std::cout << "Running " << runs[i] << " grid points" << std::endl;
             for (int reps = 0; reps < 5; reps++) {
                 timer.start(i);
-                PFHub1aPeriodic simul(runs[i], 500);
-                simul.timestep(500);
+                PFHub1aPeriodic simul(runs[i], dt);
+                simul.run_until_time(end_time);
                 timer.stop(i);
             }
         }

examples/PFHub1a.cpp

@@ -0,0 +1,34 @@
+#include <PFHub.hpp>
+
+using namespace CabanaPF;
+
+int main(int argc, char* argv[]) {
+    MPI_Init(&argc, &argv);
+    {
+        Kokkos::ScopeGuard scope_guard(argc, argv);
+        try {
+            if (argc != 5)
+                throw std::invalid_argument("");
+            const int grid_points = std::stoi(argv[2]);
+            const double dt = std::stod(argv[3]);
+            const double end_time = std::stod(argv[4]);
+
+            // read in the type of simulation to run:
+            std::unique_ptr<PFHub1aBase> simulation;
+            std::string problem_name(argv[1]);
+            if (problem_name == "benchmark")
+                simulation = std::make_unique<PFHub1aBenchmark>(grid_points, dt);
+            else if (problem_name == "periodic")
+                simulation = std::make_unique<PFHub1aPeriodic>(grid_points, dt);
+            else
+                throw std::invalid_argument("");
+
+            simulation->run_until_time(end_time);
+            simulation->output();
+        } catch (std::invalid_argument const&) {
+            std::cout << "Usage: ./PFHub1a <benchmark|periodic> <grid points> <dt> <end time>" << std::endl;
+        }
+    }
+    MPI_Finalize();
+    return 0;
+}

src/PFHub.hpp

@@ -16,24 +16,23 @@ class PFHub1aBase : public CabanaPFRunner<2> {
   protected:
     using cdouble = Kokkos::complex<double>;
 
-    const double cell_size;
-    const int timesteps;
-    const int grid_points;
     Kokkos::View<cdouble**, device_type> laplacian_view;
+    PFVariables<2, 2> vars;
 
   public:
-    PFVariables<2, 2> vars;
     static constexpr double _SIZE = 200.;
-    static constexpr double END_TIME = 250.;
     static constexpr double _KAPPA = 2.;
     static constexpr double _M = 5.;
     static constexpr double _RHO = 5.;
     static constexpr double _C_ALPHA = .3;
     static constexpr double _C_BETA = .7;
 
-    PFHub1aBase(int grid_points, int timesteps)
-        : CabanaPFRunner(grid_points, timesteps, _SIZE), cell_size{_SIZE / grid_points}, timesteps{timesteps},
-          grid_points{grid_points}, vars{layout, {"c", "df_dc"}} {
+    const int grid_points;
+    const double cell_size;
+
+    PFHub1aBase(int grid_points, double dt)
+        : CabanaPFRunner(grid_points, _SIZE, dt), vars{layout, {"c", "df_dc"}},
+          grid_points{grid_points}, cell_size{_SIZE / grid_points} {
         laplacian_view = Kokkos::View<cdouble**, device_type>("laplacian", grid_points, grid_points);
     }
 
@@ -68,7 +67,7 @@ class PFHub1aBase : public CabanaPFRunner<2> {
         initial_conditions();
     }
 
-    void pre_step() override {
+    void calc_dfdc() {
         // Calculate df_dc values:
         const auto c_view = vars[0];
         const auto dfdc_view = vars[1];
@@ -85,16 +84,16 @@ class PFHub1aBase : public CabanaPFRunner<2> {
     }
 
     void step() override {
+        calc_dfdc();
         // enter Fourier space:
         vars.fft_forward(0);
         vars.fft_forward(1);
 
-        const double dt = END_TIME / timesteps;
-        const double M = _M, KAPPA = _KAPPA;
+        // step forward with semi-implicit Euler in Fourier space:
+        const double dt = this->dt, M = _M, KAPPA = _KAPPA;
         const auto c = vars[0];
         const auto df_dc = vars[1];
         const auto laplacian = laplacian_view;
-
         node_parallel_for(
             "timestep", KOKKOS_LAMBDA(const int i, const int j) {
                 const cdouble df_dc_hat(df_dc(i, j, 0), df_dc(i, j, 1));
@@ -105,12 +104,15 @@ class PFHub1aBase : public CabanaPFRunner<2> {
                 c(i, j, 0) = c_hat.real();
                 c(i, j, 1) = c_hat.imag();
             });
-    }
-
-    void post_step() override {
-        // rescue concentration values from Fourier space:
+        // rescue concentration values from Fourier space (dfdc can be left there since it gets recalculated next
+        // timestep anyways)
         vars.fft_inverse(0);
     }
+
+    virtual void output() {}
+
+    // needed to allow polymorphism:
+    virtual ~PFHub1aBase() {}
 };
 
 class PFHub1aBenchmark : public PFHub1aBase {
@@ -131,13 +133,13 @@ class PFHub1aBenchmark : public PFHub1aBase {
             });
     }
 
-    void output() {
+    void output() override {
         std::stringstream s;
-        s << "1aBenchmark_N" << grid_points << "T" << timesteps;
+        s << "1aBenchmark_N" << grid_points << "_DT" << std::fixed << std::setprecision(3) << std::scientific << dt;
         vars.save(0, s.str());
     }
 
-    PFHub1aBenchmark(int grid_points, int timesteps) : PFHub1aBase{grid_points, timesteps} {}
+    PFHub1aBenchmark(int grid_points, double dt) : PFHub1aBase{grid_points, dt} {}
 };
 
 class PFHub1aPeriodic : public PFHub1aBase {
@@ -159,13 +161,13 @@ class PFHub1aPeriodic : public PFHub1aBase {
             });
     }
 
-    void output() {
+    void output() override {
         std::stringstream s;
-        s << "1aPeriodic_N" << grid_points << "T" << timesteps;
+        s << "1aPeriodic_N" << grid_points << "_DT" << std::fixed << std::setprecision(3) << std::scientific << dt;
         vars.save(0, s.str());
     }
 
-    PFHub1aPeriodic(int grid_points, int timesteps) : PFHub1aBase{grid_points, timesteps} {}
+    PFHub1aPeriodic(int grid_points, double dt) : PFHub1aBase{grid_points, dt} {}
 };
 
 } // namespace CabanaPF

src/Runner.hpp

@@ -20,10 +20,10 @@ class CabanaPFRunner {
 
   public:
     const int grid_points;
-    const int timesteps;
+    const double dt;
 
-    CabanaPFRunner(int grid_points, int timesteps, double size)
-        : timesteps_done{0}, have_initialized{false}, grid_points{grid_points}, timesteps{timesteps} {
+    CabanaPFRunner(int grid_points, double size, double dt)
+        : timesteps_done{0}, have_initialized{false}, grid_points{grid_points}, dt{dt} {
         std::array<double, NumSpaceDim> low_corner;
         low_corner.fill(0.0);
         std::array<double, NumSpaceDim> high_corner;
@@ -55,27 +55,39 @@ class CabanaPFRunner {
         return result;
     }
 
-    void timestep(int count) {
+    // Do a certain number of timesteps:
+    void run(const int timesteps) {
         if (!have_initialized) {
             initialize();
             have_initialized = true;
         }
-        for (int i = 0; i < count; i++) {
-            pre_step();
+        for (int i = 0; i < timesteps; i++) {
             step();
-            post_step();
             timesteps_done++;
         }
-        if (timesteps_done == timesteps)
-            finalize();
     }
 
-    // Generally, you inherit from this class and implement one or more of these:
+    void run_for_time(const double time) {
+        run(round(time / dt));
+    }
+
+    // runs as many timesteps as are needed to have done a certain number
+    void run_until_steps(const int timesteps) {
+        run(timesteps - timesteps_done);
+    }
+
+    // runs until a certain time.  Denominator allows for fractional times
+    void run_until_time(const double time) {
+        run(round(time / dt - timesteps_done));
+    }
+
+    int get_timesteps_done() const {
+        return timesteps_done;
+    }
+
+    // Children inherit from this class and implement these:
     virtual void initialize() {} // Called once, before taking the first timestep
-    virtual void pre_step() {}   // Called each timestep
-    virtual void step() {}       // Called each timestep, after pre_step
-    virtual void post_step() {}  // Called each timestep, after step
-    virtual void finalize() {}   // Called once, when the requested number of timesteps have been done
+    virtual void step() {}       // Called to take a timestep
 };
 
 } // namespace CabanaPF

unit_test/test.cpp

@@ -14,8 +14,8 @@ PFHub1a Available here (ORNL internal): https://code.ornl.gov/71d/phase-field-ex
 selected
 */
 TEST(PFHub1a, Initialization) {
-    PFHub1aBenchmark simulation(96, 500);
-    simulation.timestep(0); // trigger initialization
+    PFHub1aBenchmark simulation(96, .5);
+    simulation.run(0); // trigger initialization
     auto results = simulation.get_cpu_view();
     //"true results" come from python implementation (see previous comment)
     // check 4 points for basic indexing/results:
@@ -38,8 +38,8 @@ TEST(PFHub1a, Initialization) {
 }
 
 TEST(PFHub1a, OneTimestep) {
-    PFHub1aBenchmark simulation(96, 500);
-    simulation.timestep(1);
+    PFHub1aBenchmark simulation(96, .5);
+    simulation.run(1);
     auto results = simulation.get_cpu_view();
     // test at extreme points and 10 random points.  Correct values come from python implementation (see above)
     EXPECT_DOUBLE_EQ(0.5214689225639189, results(0, 0, 0));
@@ -57,8 +57,8 @@ TEST(PFHub1a, OneTimestep) {
 }
 
 TEST(PFHub1a, AllTimestep) {
-    PFHub1aBenchmark simulation(96, 500);
-    simulation.timestep(500);
+    PFHub1aBenchmark simulation(96, .5);
+    simulation.run(500);
     auto results = simulation.get_cpu_view();
     // as before, (0,0), (95,95), and 10 random points, testing against python
     EXPECT_NEAR(0.4412261765305555, results(0, 0, 0), 1e-9);
@@ -106,24 +106,24 @@ TEST(PFVariables, saveload) {
 
 // Similar to above, the python implementation was modified to use the periodic initial conditions
 TEST(PFHub1aPeriodic, periodic) {
-    PFHub1aPeriodic simulation(96, 500);
-    simulation.timestep(0);
+    PFHub1aPeriodic simulation(96, .5);
+    simulation.run(0);
     auto results = simulation.get_cpu_view();
     EXPECT_NEAR(0.53, results(0, 0, 0), 1e-9);
     EXPECT_NEAR(0.5280872555770657, results(95, 95, 0), 1e-9);
     EXPECT_NEAR(0.49625, results(56, 52, 0), 1e-9);
     EXPECT_NEAR(0.5096103712433676, results(39, 36, 0), 1e-9);
     EXPECT_NEAR(0.5122826024701564, results(46, 40, 0), 1e-9);
 
-    simulation.timestep(1);
+    simulation.run_until_time(.5);
     results = simulation.get_cpu_view();
     EXPECT_NEAR(0.5316722086053631, results(0, 0, 0), 1e-9);
     EXPECT_NEAR(0.5296339912527902, results(95, 95, 0), 1e-9);
     EXPECT_NEAR(0.5155424558547776, results(24, 46, 0), 1e-9);
     EXPECT_NEAR(0.510243048825588, results(87, 78, 0), 1e-9);
     EXPECT_NEAR(0.5092351158827323, results(6, 19, 0), 1e-9);
 
-    simulation.timestep(499);
+    simulation.run_until_steps(500);
     results = simulation.get_cpu_view();
     EXPECT_NEAR(0.6993369106233298, results(0, 0, 0), 1e-9);
     EXPECT_NEAR(0.7014658707445363, results(95, 95, 0), 1e-9);