Low storage runge kutta #259

jocelynpoon · 2024-06-17T15:00:50Z

Draft PR for the low storage RK method code

sivanadarajah · 2024-06-21T17:15:30Z

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

+    //int k = q_hat +1;
+
+    double sum_delta = 0;
+    double atol = 0.1;


Lets include atol and rtol as parameters

sivanadarajah · 2024-06-21T17:16:07Z

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

+    int k = q_hat +1;
+   // double error = 0.0;
+    //double w = 0.0;
+    double beta_controller[3] = {0.70, -0.40, 0};


The same for beta as well.

cpethrick

This is great work! I'm glad to see it coming together so well. I've left some preliminary comments to work through while the PR is still in draft. As you clean up your code, it is also good to remove commented lines, though I know those are sometimes necessary while you're still developing it!

cpethrick · 2024-07-12T15:53:28Z

src/flow_solver/flow_solver.cpp


            // Compute the unsteady quantities, write to the dealii table, and output to file
            flow_solver_case->compute_unsteady_data_and_write_to_table(ode_solver, dg, unsteady_data_table);
            // update next time step
-            if(flow_solver_param.adaptive_time_step == true) {


Could you add another parameter which controls whether you use the error-adaptive time stepping? Then this structure would change to

if (CFL adaptive timestep) next_time_step = flow_solver_case->get_adaptive_time_step(dg) else if (error adaptive time step) next_time_step = ode_solver->err_time_step(time_step,false); else next_time_step = flow_solver_case->get_constant_time_step(dg);

I think the best place for that parameter would be in the flow solver section.

cpethrick · 2024-07-12T15:55:03Z

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

+}
+
+template <int dim, typename real, int n_rk_stages, typename MeshType>
+void LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::step_in_time (real dt, const bool pseudotime)


Suggested change

void LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::step_in_time (real dt, const bool pseudotime)

void LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::step_in_time (real dt, const bool /*pseudotime*/)

This way you don't have to call (void) pseudotime.

cpethrick · 2024-07-12T15:59:04Z

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

+    dealii::LinearAlgebra::distributed::Vector<double> storage_register_1;
+    storage_register_1.reinit(this->solution_update);
+    storage_register_1 = this->dg->solution;
+    dealii::LinearAlgebra::distributed::Vector<double> storage_register_3 = storage_register_1;


Could you move the initialization of these vectors to allocate_ode_system()?

In order to do that, you would need to make all of the storage registers member variables of the LowStorageRungeKuttaODESolver, then just do the allocation steps in allocate_ode_system().

Just be careful to reinitialize the vectors if that needs to be done!

cpethrick · 2024-07-12T16:00:48Z

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

+    for (int i = 1; i < n_rk_stages +1; i++ ){
+        // storage_register_2 = storage_register_2 + delta[i-1] * storage_register_1;
+        storage_register_2.add(this->butcher_tableau->get_delta(i-1) , storage_register_1);
+        this->dg->solution = rhs;


I think this might be possible to do without the extra vector rhs, and instead using storage_register_1. I haven't gone through the loop in detail, but maybe you could check that.

cpethrick · 2024-07-12T16:02:50Z

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

+    }
+    this->dg->solution = storage_register_1;
+
+    double global_size = dealii::Utilities::MPI::sum(storage_register_1.local_size(), this->mpi_communicator);


This would probably be better to do in allocate_ode_system and storing as a member variable. Summing over all of the cores is often quite slow, especially as you increase the number of parallel processes.

cpethrick · 2024-07-12T16:20:59Z

src/ode_solver/low_storage_runge_kutta_ode_solver.h

+    /** This is initialized for any RK method, but solution-sized vectors are 
+     *  only initialized if there is an implicit solve
+     */
+    JFNKSolver<dim,real,MeshType> solver;


I think this can be removed since we will never have implicit for this type of low-storage algorithm. You can remove it from the header here and also from line 13 of the cpp file.

cpethrick · 2024-07-12T16:24:15Z

src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.h

+namespace PHiLiP {
+namespace ODE {
+
+#if PHILIP_DIM==1


Could you also add references to the paper and github repo to this class?

cpethrick · 2024-07-12T16:25:15Z

src/parameters/parameters_linear_solver.cpp

@@ -7,6 +7,18 @@ void LinearSolverParam::declare_parameters (dealii::ParameterHandler &prm)
 {
    prm.enter_subsection("linear solver");
    {
+        prm.declare_entry("atol", "0.001",


Could these be moved to the ode_solver parameters section rather than this one? I think they would be more relevant there.

cpethrick · 2024-07-12T16:26:10Z

src/parameters/parameters_ode_solver.cpp

-                          " dirk_3_im"),
+                          " dirk_3_im | "
+                          " RK3_2_5F_3SStarPlus | "
+                          " RK4_3_5_3SStar "),
                          "Runge-kutta method to use. Methods with _ex are explicit, and with _im are implicit."


Maybe add another line to this description indicating that 2SStarPlus and 3SStar are low-storage RK methods.

cpethrick · 2024-07-12T16:27:03Z

src/parameters/parameters_ode_solver.h

@@ -61,12 +62,15 @@ class ODESolverParam
        euler_ex, ///Forward Euler
        euler_im, ///Implicit Euler
        dirk_2_im, ///Second-order diagonally-implicit RK
-        dirk_3_im ///Third-order diagonally-implicit RK
+        dirk_3_im, ///Third-order diagonally-implicit RK
+        RK4_3_5_3SStar,


Add short comments with the main properties - include order, number of stages, FSAL, and other pertinent information.

jocelynpoon added 9 commits June 12, 2024 11:37

First commit of new Low Storage RK ODESolver. [Does not compile]

4f19280

low storage test

93f71f1

compiles

01d72b3

this works

ca60d55

compiles

555565c

compiles but time step does not change as expected

8cb35d8

fixed w term need to fix dt

303cbf6

this runs

d858fe7

step size change

35170cc

sivanadarajah reviewed Jun 21, 2024

View reviewed changes

jocelynpoon added 9 commits June 25, 2024 10:33

save

15ad87c

works with mpi > 1

820d80c

update

e758b4a

add test euler inviscid tgv largest stable cfl

ee073a4

add 3s*+ method and RK3(2)5F[3S*+]

7acd254

low storage algorithms, need to fix coefficients for RK3_2_5F_3SStarPlus

40d4e00

correct convergence for both algorithms (3S* and 3S*+)

659719b

add tolerances as parameters

4a99eac

tolerance parameters work

ac2b446

cpethrick reviewed Jul 12, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Low storage runge kutta #259

Low storage runge kutta #259

jocelynpoon commented Jun 17, 2024

sivanadarajah Jun 21, 2024

sivanadarajah Jun 21, 2024

cpethrick left a comment

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

cpethrick Jul 12, 2024

	void LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::step_in_time (real dt, const bool pseudotime)
	void LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::step_in_time (real dt, const bool /pseudotime/)

Low storage runge kutta #259

Are you sure you want to change the base?

Low storage runge kutta #259

Conversation

jocelynpoon commented Jun 17, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cpethrick left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment