Framework for adaptive dt / independent dt specs in each solver input (#107)

* framework for adaptive dt / independent dt specs in each solver input

* put adaptive dt into a reg test

* change function call to go with amr-wind changes

* get min dt across processors

* check for fixed dt upon initialization; do other work for consistent dt only when needed

* allow max_time as ending argument, modify reg test to include it

* remove redundant mpi_barrier

---------

Co-authored-by: Jon Rood <[email protected]>

Author: mbkuhn

app/exawind/exawind.cpp

@@ -186,7 +186,10 @@ int main(int argc, char** argv)
     sim.set_nw_start_rank(nalu_start_rank);
 
     const auto nalu_vars = node["nalu_vars"].as<std::vector<std::string>>();
-    const int num_timesteps = node["num_timesteps"].as<int>();
+    const int num_timesteps =
+        node["num_timesteps"] ? node["num_timesteps"].as<int>() : -1;
+    const double max_time =
+        node["max_time"] ? node["max_time"].as<double>() : -1.0;
     const int additional_picard_its =
         node["additional_picard_iterations"]
             ? node["additional_picard_iterations"].as<int>()
@@ -199,6 +202,12 @@ int main(int argc, char** argv)
                                  : false;
     sim.set_holemap_alg(holemap_alg);
 
+    if (num_timesteps < 0 && max_time < 0.) {
+        throw std::runtime_error(
+            "max_timesteps or num_timesteps must be specified as positive "
+            "values. These are both unspecified or specified as negative.");
+    }
+
     if (node["composite_body"]) {
         const YAML::Node& composite_mesh = node["composite_body"];
         const int num_composite = static_cast<int>(composite_mesh.size());
@@ -291,7 +300,8 @@ int main(int argc, char** argv)
     sim.echo("Initializing overset simulation");
     sim.initialize();
     sim.echo("Initialization successful");
-    sim.run_timesteps(additional_picard_its, nonlinear_its, num_timesteps);
+    sim.run_timesteps(
+        additional_picard_its, nonlinear_its, num_timesteps, max_time);
     sim.delete_solvers();
 
     if (amr_comm != MPI_COMM_NULL) {

src/AMRWind.cpp

@@ -70,11 +70,18 @@ void AMRWind::prepare_solver_epilog()
     m_incflo.prepare_for_time_integration();
 }
 
-void AMRWind::pre_advance_stage1(size_t inonlin)
+void AMRWind::pre_advance_stage0(size_t inonlin)
 {
     if (inonlin < 1) {
         m_incflo.sim().time().new_timestep();
         m_incflo.regrid_and_update();
+        m_incflo.compute_dt();
+    }
+}
+
+void AMRWind::pre_advance_stage1(size_t inonlin)
+{
+    if (inonlin < 1) {
         m_incflo.pre_advance_stage1();
     }
 }
@@ -84,6 +91,20 @@ void AMRWind::pre_advance_stage2(size_t inonlin)
     if (inonlin < 1) m_incflo.pre_advance_stage2();
 }
 
+double AMRWind::get_time() { return m_incflo.time().new_time(); }
+
+double AMRWind::get_timestep_size() { return m_incflo.time().delta_t(); }
+
+void AMRWind::set_timestep_size(const double dt)
+{
+    m_incflo.sim().time().delta_t() = dt;
+}
+
+bool AMRWind::is_fixed_timestep_size()
+{
+    return (!m_incflo.sim().time().adaptive_timestep());
+}
+
 void AMRWind::advance_timestep(size_t inonlin) { m_incflo.do_advance(inonlin); }
 
 void AMRWind::post_advance() { m_incflo.post_advance_work(); }

src/AMRWind.h

@@ -27,6 +27,7 @@ class AMRWind : public ExawindSolver
     ~AMRWind();
     bool is_unstructured() override { return false; }
     bool is_amr() override { return true; }
+    bool is_fixed_timestep_size() override;
     int overset_update_interval() override;
     int time_index() override;
     std::string identifier() override { return "AMR-Wind"; }
@@ -37,8 +38,12 @@ class AMRWind : public ExawindSolver
     void init_epilog() override;
     void prepare_solver_prolog() override;
     void prepare_solver_epilog() override;
+    void pre_advance_stage0(size_t inonlin) override;
     void pre_advance_stage1(size_t inonlin) override;
     void pre_advance_stage2(size_t inonlin) override;
+    double get_time() override;
+    double get_timestep_size() override;
+    void set_timestep_size(const double) override;
     void advance_timestep(size_t inonlin) override;
     void additional_picard_iterations(const int) override {};
     void post_advance() override;

src/ExawindSolver.h

@@ -19,6 +19,13 @@ class ExawindSolver
     void call_init_epilog() { init_epilog(); };
     void call_prepare_solver_prolog() { prepare_solver_prolog(); };
     void call_prepare_solver_epilog() { prepare_solver_epilog(); };
+    void call_pre_advance_stage0(size_t inonlin, const bool increment)
+    {
+        const std::string name = "Pre";
+        m_timers.tick(name, increment);
+        pre_advance_stage0(inonlin);
+        m_timers.tock(name);
+    }
     void call_pre_advance_stage1(size_t inonlin, const bool increment)
     {
         const std::string name = "Pre";
@@ -33,6 +40,29 @@ class ExawindSolver
         pre_advance_stage2(inonlin);
         m_timers.tock(name);
     };
+    double call_get_time()
+    {
+        const std::string name = "Pre";
+        m_timers.tick(name);
+        double time = get_time();
+        m_timers.tock(name);
+        return time;
+    }
+    double call_get_timestep_size()
+    {
+        const std::string name = "Pre";
+        m_timers.tick(name);
+        double dt = get_timestep_size();
+        m_timers.tock(name);
+        return dt;
+    };
+    void call_set_timestep_size(double dt)
+    {
+        const std::string name = "Pre";
+        m_timers.tick(name);
+        set_timestep_size(dt);
+        m_timers.tock(name);
+    };
     void call_advance_timestep(size_t inonlin, const bool increment)
     {
         const std::string name = "Solve";
@@ -91,6 +121,7 @@ class ExawindSolver
 
     virtual bool is_unstructured() { return false; };
     virtual bool is_amr() { return false; };
+    virtual bool is_fixed_timestep_size() = 0;
     virtual int overset_update_interval() { return 100000000; };
     virtual int time_index() = 0;
     virtual std::string identifier() { return "ExawindSolver"; }
@@ -108,8 +139,12 @@ class ExawindSolver
     virtual void init_epilog() = 0;
     virtual void prepare_solver_prolog() = 0;
     virtual void prepare_solver_epilog() = 0;
+    virtual void pre_advance_stage0(size_t inonlin) = 0;
     virtual void pre_advance_stage1(size_t inonlin) = 0;
     virtual void pre_advance_stage2(size_t inonlin) = 0;
+    virtual double get_time() = 0;
+    virtual double get_timestep_size() = 0;
+    virtual void set_timestep_size(const double) = 0;
     virtual void advance_timestep(size_t inonlin) = 0;
     virtual void additional_picard_iterations(const int) = 0;
     virtual void post_advance() = 0;

src/NaluWind.cpp

@@ -82,6 +82,11 @@ void NaluWind::prepare_solver_epilog()
         realm->output_converged_results();
 }
 
+void NaluWind::pre_advance_stage0(size_t inonlin)
+{
+    m_sim.timeIntegrator_->prepare_time_step(inonlin);
+}
+
 void NaluWind::pre_advance_stage1(size_t inonlin)
 {
     m_sim.timeIntegrator_->pre_realm_advance_stage1(inonlin);
@@ -92,6 +97,26 @@ void NaluWind::pre_advance_stage2(size_t inonlin)
     m_sim.timeIntegrator_->pre_realm_advance_stage2(inonlin);
 }
 
+double NaluWind::get_time()
+{
+    return m_sim.timeIntegrator_->get_current_time();
+}
+
+double NaluWind::get_timestep_size()
+{
+    return m_sim.timeIntegrator_->get_time_step();
+}
+
+void NaluWind::set_timestep_size(const double dt)
+{
+    m_sim.timeIntegrator_->set_time_step(dt);
+}
+
+bool NaluWind::is_fixed_timestep_size()
+{
+    return m_sim.timeIntegrator_->get_is_fixed_time_step();
+}
+
 void NaluWind::advance_timestep(size_t /*inonlin*/)
 {
     for (auto* realm : m_sim.timeIntegrator_->realmVec_) {

src/NaluWind.h

@@ -50,6 +50,7 @@ class NaluWind : public ExawindSolver
     ~NaluWind();
     bool is_unstructured() override { return true; }
     bool is_amr() override { return false; }
+    bool is_fixed_timestep_size() override;
     int overset_update_interval() override;
     int time_index() override;
     std::string identifier() override
@@ -64,8 +65,12 @@ class NaluWind : public ExawindSolver
     void init_epilog() override;
     void prepare_solver_prolog() override;
     void prepare_solver_epilog() override;
+    void pre_advance_stage0(size_t inonlin) override;
     void pre_advance_stage1(size_t inonlin) override;
     void pre_advance_stage2(size_t inonlin) override;
+    double get_time() override;
+    double get_timestep_size() override;
+    void set_timestep_size(const double) override;
     void advance_timestep(size_t inonlin) override;
     void additional_picard_iterations(const int) override;
     void post_advance() override;

src/OversetSimulation.cpp

@@ -82,7 +82,12 @@ void OversetSimulation::initialize()
     }
     m_last_timestep = m_solvers.at(0)->time_index();
 
-    MPI_Barrier(m_comm);
+    // Determine if any of the solvers have adaptive timestepping
+    for (auto& ss : m_solvers) {
+        m_fixed_dt = m_fixed_dt && ss->is_fixed_timestep_size();
+    }
+    MPI_Allreduce(MPI_IN_PLACE, &m_fixed_dt, 1, MPI_C_BOOL, MPI_LAND, m_comm);
+
     m_initialized = true;
 }
 
@@ -127,7 +132,10 @@ void OversetSimulation::exchange_solution(bool increment_time)
 }
 
 void OversetSimulation::run_timesteps(
-    const int add_pic_its, const int nonlinear_its, const int nsteps)
+    const int add_pic_its,
+    const int nonlinear_its,
+    const int nsteps,
+    const double max_time)
 {
 
     if (!m_initialized) {
@@ -140,7 +148,12 @@ void OversetSimulation::run_timesteps(
         "Running " + std::to_string(nsteps) + " timesteps starting from " +
         std::to_string(tstart));
 
-    for (int nt = tstart; nt < tend; ++nt) {
+    int nt = tstart;
+    double time = max_time < 0. ? 0. : m_solvers[0]->call_get_time();
+    bool step_check = nsteps > 0 ? nt < tend : true;
+    bool time_check = max_time > 0. ? time < max_time : true;
+    bool do_step = step_check && time_check;
+    while (do_step) {
         m_printer.echo_time_header();
 
         m_timers_exa.tick("TimeStep");
@@ -150,6 +163,20 @@ void OversetSimulation::run_timesteps(
 
             bool increment_timer = inonlin > 0 ? true : false;
 
+            double dt{1e8};
+            for (auto& ss : m_solvers) {
+                ss->call_pre_advance_stage0(inonlin, increment_timer);
+                if (inonlin < 1 && !m_fixed_dt) {
+                    dt = std::min(dt, ss->call_get_timestep_size());
+                }
+            }
+
+            if (inonlin < 1 && !m_fixed_dt) {
+                MPI_Allreduce(
+                    MPI_IN_PLACE, &dt, 1, MPI_DOUBLE, MPI_MIN, m_comm);
+                for (auto& ss : m_solvers) ss->call_set_timestep_size(dt);
+            }
+
             for (auto& ss : m_solvers)
                 ss->call_pre_advance_stage1(inonlin, increment_timer);
 
@@ -183,6 +210,12 @@ void OversetSimulation::run_timesteps(
         MPI_Barrier(m_comm);
 
         mem_usage_all(nt);
+
+        ++nt;
+        if (max_time > 0.) time = m_solvers[0]->call_get_time();
+        step_check = nsteps > 0 ? nt < tend : true;
+        time_check = max_time > 0. ? time < max_time : true;
+        do_step = step_check && time_check;
     }
     for (auto& ss : m_solvers) ss->call_dump_simulation_time();
     m_last_timestep = tend;

src/OversetSimulation.h

@@ -30,6 +30,9 @@ class OversetSimulation
     int m_overset_update_interval{100000000};
     //! Last timestep run during this simulation
     int m_last_timestep{0};
+    //! Flag indicating whether all solvers use fixed dt. If any solver uses
+    //! adaptive dt, then this flag will be false
+    bool m_fixed_dt{true};
     //! Flag indicating whether initialization tasks have been performed
     bool m_initialized{false};
     //! Flag indicating if complementary comms have been initialized
@@ -88,7 +91,8 @@ class OversetSimulation
     void run_timesteps(
         const int add_pic_its,
         const int nonlinear_its = 1,
-        const int nsteps = 1);
+        const int nsteps = 1,
+        const double max_time = -1.);
 
     //! Print something
     void echo(const std::string& out) { m_printer.echo(out); }

test/test_files/dam-break-block/dam-break-block-amr.inp

@@ -1,8 +1,8 @@
 time.stop_time               =   2     # Max (simulated) time to evolve
 time.max_step                =   2000          # Max number of time steps
 
-time.fixed_dt         =   0.005        # Use this constant dt if > 0
-time.cfl              =   1.0         # CFL factor
+time.initial_dt       =   0.01        # Use this constant dt if > 0
+time.cfl              =   0.95         # CFL factor
 time.use_force_cfl    =   false
 
 time.plot_interval            =  10       # Steps between plot files

test/test_files/dam-break-block/dam-break-block-nalu.yaml

@@ -33,6 +33,10 @@ realms:
     automatic_decomposition_type: rcb
     check_for_missing_bcs: yes
 
+    time_step_control:
+     target_courant: 2.0
+     time_step_change_factor: 1.2
+
     equation_systems:
       name: theEqSys
       max_iterations: 3
@@ -132,8 +136,8 @@ Time_Integrators:
       name: ti_1
       start_time: 0.0
       termination_time: 2.0
-      time_step: 0.005
-      time_stepping_type: fixed
+      time_step: 0.01
+      time_stepping_type: adaptive
       time_step_count: 0
       second_order_accuracy: yes