Low storage runge kutta

.gitignore

@@ -26,3 +26,4 @@ tests/integration_tests_control_files/decaying_homogeneous_isotropic_turbulence/
 tests/integration_tests_control_files/decaying_homogeneous_isotropic_turbulence/setup_files/4proc/*.dat
 tests/unit_tests/grid/gmsh_reader/*.msh
 tests/meshes/*.msh
+*.vscode

src/flow_solver/flow_solver.cpp

@@ -500,19 +500,23 @@ int FlowSolver<dim,nstate>::run() const
             }
 
             // update time step in flow_solver_case
-            flow_solver_case->set_time_step(time_step);
-
+            ode_solver->step_in_time(time_step,false);
+            pcout << "time step " << time_step;
+            next_time_step = ode_solver->err_time_step(time_step,false); 
+            // pcout << "next time step " << next_time_step;
             // advance solution
-            ode_solver->step_in_time(time_step,false); // pseudotime==false
+            // ode_solver->err_time_step(time_step); // pseudotime==false
 
             // 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) {
+            /*if(flow_solver_param.adaptive_time_step == true) {
                 next_time_step = flow_solver_case->get_adaptive_time_step(dg);
             } else {
                 next_time_step = flow_solver_case->get_constant_time_step(dg);
             }
+            */
+
 
 #if PHILIP_DIM>1
             if(flow_solver_param.output_restart_files == true) {

src/flow_solver/flow_solver_cases/periodic_turbulence.cpp

@@ -1,5 +1,5 @@
 #include "periodic_turbulence.h"
-
+#include "ode_solver/low_storage_runge_kutta_ode_solver.h"
 #include <deal.II/base/function.h>
 #include <stdlib.h>
 #include <iostream>
@@ -686,7 +686,15 @@ void PeriodicTurbulence<dim, nstate>::compute_unsteady_data_and_write_to_table(
     //unpack current iteration and current time from ode solver
     const unsigned int current_iteration = ode_solver->current_iteration;
     const double current_time = ode_solver->current_time;
-
+    //std::shared_ptr<ODE::LowStorageRungeKuttaODESolver<dim, double, 5>> ode_solver_low_storage = std::dynamic_pointer_cast<ODE::LowStorageRungeKuttaODESolver<dim, double, 5>>(ode_solver);
+    //(void) ode_solver_low_storage;
+    //double eps0 = ode_solver_low_storage->epsilon[0];
+    //double eps1 = ode_solver_low_storage->epsilon[1];
+    //double eps2 = ode_solver_low_storage->epsilon[2];
+    //this->pcout << eps0 << std::endl;
+    //std::abort();
+    //double eps1[3] = ode_solver_low_storage.epsilon[1];
+    //double eps2[3] = ode_solver_low_storage.epsilon[2];
     // Compute and update integrated quantities
     this->compute_and_update_integrated_quantities(*dg);
     // Get computed quantities
@@ -707,6 +715,9 @@ void PeriodicTurbulence<dim, nstate>::compute_unsteady_data_and_write_to_table(
 
     if(this->mpi_rank==0) {
         // Add values to data table
+        //this->add_value_to_data_table(eps0,"epsilon0",unsteady_data_table);
+        //this->add_value_to_data_table(eps1,"epsilon1",unsteady_data_table);
+        //this->add_value_to_data_table(eps2,"epsilon2",unsteady_data_table);
         this->add_value_to_data_table(current_time,"time",unsteady_data_table);
         if(do_calculate_numerical_entropy) this->add_value_to_data_table(this->cumulative_numerical_entropy_change_FRcorrected,"numerical_entropy_scaled_cumulative",unsteady_data_table);
         if(is_rrk) this->add_value_to_data_table(relaxation_parameter, "relaxation_parameter",unsteady_data_table);

src/ode_solver/CMakeLists.txt

@@ -2,8 +2,11 @@ set(ODE_SOURCE
     ode_solver_factory.cpp
     ode_solver_base.cpp
     runge_kutta_ode_solver.cpp
+    low_storage_runge_kutta_ode_solver.cpp
     runge_kutta_methods/runge_kutta_methods.cpp
     runge_kutta_methods/rk_tableau_base.cpp
+    runge_kutta_methods/low_storage_rk_tableau_base.cpp
+    runge_kutta_methods/low_storage_runge_kutta_methods.cpp
     relaxation_runge_kutta/empty_RRK_base.cpp
     relaxation_runge_kutta/runge_kutta_store_entropy.cpp
     relaxation_runge_kutta/rrk_ode_solver_base.cpp

src/ode_solver/low_storage_runge_kutta_ode_solver.cpp

@@ -0,0 +1,188 @@
+#include "low_storage_runge_kutta_ode_solver.h"
+
+namespace PHiLiP {
+namespace ODE {
+
+template <int dim, typename real, int n_rk_stages, typename MeshType> 
+LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::LowStorageRungeKuttaODESolver(std::shared_ptr< DGBase<dim, real, MeshType> > dg_input,
+        std::shared_ptr<LowStorageRKTableauBase<dim,real,MeshType>> rk_tableau_input,
+        std::shared_ptr<EmptyRRKBase<dim,real,MeshType>> RRK_object_input)
+        : ODESolverBase<dim,real,MeshType>(dg_input)
+        , butcher_tableau(rk_tableau_input)
+        , relaxation_runge_kutta(RRK_object_input)
+        , solver(dg_input)
+{       epsilon[0] = 1.0;
+        epsilon[1] = 1.0;
+        epsilon[2] = 1.0; 
+}
+
+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)
+{
+    this->original_time_step = dt;
+    this->solution_update = this->dg->solution; //storing u_n
+    (void) pseudotime;
+
+    dealii::LinearAlgebra::distributed::Vector<double> storage_register_2;
+    storage_register_2.reinit(this->solution_update);
+    storage_register_2*=0;
+    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;
+
+    dealii::LinearAlgebra::distributed::Vector<double> rhs = storage_register_1;
+
+    dealii::LinearAlgebra::distributed::Vector<double> storage_register_4;
+
+    if (this->butcher_tableau->get_b_hat(0) != 0.0){
+        storage_register_4 = storage_register_1;
+    } 
+
+    double sum_delta = 0;
+
+    //double atol = 0.001;
+    //double rtol = 0.001;
+    double error = 0.0;
+    w = 0.0;
+
+    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;
+        this->dg->assemble_residual();
+        this->dg->apply_inverse_global_mass_matrix(this->dg->right_hand_side, rhs);
+        // storage_register_1 = gamma[i][0] * storage_register_1 + gamma[i][1] * storage_register_2 + gamma[i][2] * storage_register_3 + beta[i] * dt * rhs;
+        storage_register_1 *= this->butcher_tableau->get_gamma(i, 0);
+        storage_register_1.add(this->butcher_tableau->get_gamma(i, 1), storage_register_2);
+        storage_register_1.add(this->butcher_tableau->get_gamma(i, 2), storage_register_3);
+        rhs *= dt;
+        storage_register_1.add(this->butcher_tableau->get_beta(i), rhs);
+        if (this->butcher_tableau->get_b_hat(i) != 0.0){
+            storage_register_4.add(this->butcher_tableau->get_b_hat(i-1), rhs);
+        }
+        //this->pcout << std::endl;
+        // Check  bhat (i)
+        // rhs = dt * f(S1)
+        rhs = storage_register_1;
+    }
+   // std::abort();
+    // storage_register_2 = (storage_register_2 + delta[m] * storage_register_1 + delta[m+1] * storage_register_3) / sum_delta;
+    if (this->butcher_tableau->get_b_hat(1) == 0){
+        for (int i = 0; i < num_delta; i++){ //change n_rk_stages+2 to num_delta on this line
+        sum_delta = sum_delta + this->butcher_tableau->get_delta(i);
+        }
+        storage_register_2.add(this->butcher_tableau->get_delta(n_rk_stages), storage_register_1);
+        storage_register_2.add(this->butcher_tableau->get_delta(n_rk_stages+1), storage_register_3);
+        storage_register_2 /= sum_delta;
+        // u_hat = s2
+    }
+
+
+    // need to calculate rhs of s1
+    else if (this->butcher_tableau->get_b_hat(n_rk_stages) != 0){
+        this->dg->solution = rhs;
+        this->dg->assemble_residual();
+        this->dg->apply_inverse_global_mass_matrix(this->dg->right_hand_side, rhs);
+        rhs *= dt;
+        storage_register_4.add(this->butcher_tableau->get_b_hat(n_rk_stages), rhs);
+        //this->pcout << " b_hat_fsal " << this->butcher_tableau->get_b_hat(n_rk_stages);
+
+    }
+    this->dg->solution = storage_register_1;
+
+    double global_size = dealii::Utilities::MPI::sum(storage_register_1.local_size(), this->mpi_communicator);
+
+    if (this->butcher_tableau->get_b_hat(1) == 0){
+        // loop sums elements at each mpi processor
+        for (dealii::LinearAlgebra::distributed::Vector<double>::size_type i = 0; i < storage_register_1.local_size(); ++i) {
+            error = storage_register_1.local_element(i) - storage_register_2.local_element(i);
+            w = w + pow(error / (atol + rtol * std::max(std::abs(storage_register_1.local_element(i)), std::abs(storage_register_2.local_element(i)))), 2);
+        }
+    }
+    else if (this->butcher_tableau->get_b_hat(1) != 0){
+        // loop sums elements at each mpi processor
+        for (dealii::LinearAlgebra::distributed::Vector<double>::size_type i = 0; i < storage_register_1.local_size(); ++i) {
+            error = storage_register_1.local_element(i) - storage_register_4.local_element(i);
+            w = w + pow(error / (atol + rtol * std::max(std::abs(storage_register_1.local_element(i)), std::abs(storage_register_4.local_element(i)))), 2);
+        }
+    }
+    this->pcout << std::endl;
+
+    // sum over all elements
+    w = dealii::Utilities::MPI::sum(w, this->mpi_communicator);
+
+    w = pow(w / global_size, 0.5);
+
+    this->pcout << std::endl;
+    // std::cout << "w " << w;
+
+    ++(this->current_iteration);
+    this->current_time += dt;
+    this->pcout << " Time is: " << this->current_time <<std::endl;
+    this->pcout << "dt" << dt;
+}
+
+
+template <int dim, typename real, int n_rk_stages, typename MeshType> 
+double LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::err_time_step (real dt, const bool pseudotime)
+{
+    (void) pseudotime;
+    //int q_hat = 2;
+    //int k = q_hat +1;
+    double beta_controller[3] = {0.70, -0.23, 0};
+
+    // error based step size 
+
+    epsilon[2] = epsilon[1];
+    epsilon[1] = epsilon[0];
+    epsilon[0] = 1.0 / w;
+    //dt = pow(epsilon[0], 1.0 * beta_controller[0]/k) * pow(epsilon[1], 1.0 * beta_controller[1]/k) * pow(epsilon[2], 1.0 * beta_controller[2]/k) * dt;
+    dt = pow(epsilon[0], 1.0 * beta_controller[0]/rk_order) * pow(epsilon[1], 1.0 * beta_controller[1]/rk_order) * pow(epsilon[2], 1.0 * beta_controller[2]/rk_order) * dt;
+    this->pcout << std::endl;
+    //this->pcout << "dt" << dt;
+    return dt;
+}
+
+
+template <int dim, typename real, int n_rk_stages, typename MeshType> 
+void LowStorageRungeKuttaODESolver<dim,real,n_rk_stages, MeshType>::allocate_ode_system ()
+{
+    this->pcout << "Allocating ODE system..." << std::flush;
+    this->solution_update.reinit(this->dg->right_hand_side);
+    if(this->all_parameters->use_inverse_mass_on_the_fly == false) {
+        this->pcout << " evaluating inverse mass matrix..." << std::flush;
+        this->dg->evaluate_mass_matrices(true); // creates and stores global inverse mass matrix
+        //RRK needs both mass matrix and inverse mass matrix
+        using ODEEnum = Parameters::ODESolverParam::ODESolverEnum;
+        ODEEnum ode_type = this->ode_param.ode_solver_type;
+        if (ode_type == ODEEnum::rrk_explicit_solver){
+            this->dg->evaluate_mass_matrices(false); // creates and stores global mass matrix
+        }
+    }
+
+    this->pcout << std::endl;
+
+    this->butcher_tableau->set_tableau();
+}
+
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,1, dealii::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,2, dealii::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,3, dealii::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,4, dealii::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,5, dealii::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,1, dealii::parallel::shared::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,2, dealii::parallel::shared::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,3, dealii::parallel::shared::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,4, dealii::parallel::shared::Triangulation<PHILIP_DIM> >;
+template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,5, dealii::parallel::shared::Triangulation<PHILIP_DIM> >;
+#if PHILIP_DIM != 1
+    template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,1, dealii::parallel::distributed::Triangulation<PHILIP_DIM> >;
+    template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,2, dealii::parallel::distributed::Triangulation<PHILIP_DIM> >;
+    template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,3, dealii::parallel::distributed::Triangulation<PHILIP_DIM> >;
+    template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,4, dealii::parallel::distributed::Triangulation<PHILIP_DIM> >;
+    template class LowStorageRungeKuttaODESolver<PHILIP_DIM, double,5, dealii::parallel::distributed::Triangulation<PHILIP_DIM> >;
+#endif
+
+} // ODESolver namespace
+} // PHiLiP namespace

src/ode_solver/low_storage_runge_kutta_ode_solver.h

@@ -0,0 +1,75 @@
+#ifndef __LOW_STORAGE_RUNGE_KUTTA_ODESOLVER__
+#define __LOW_STORAGE_RUNGE_KUTTA_ODESOLVER__
+
+#include "JFNK_solver/JFNK_solver.h"
+#include "dg/dg_base.hpp"
+#include "ode_solver_base.h"
+#include "runge_kutta_methods/low_storage_rk_tableau_base.h"
+#include "relaxation_runge_kutta/empty_RRK_base.h"
+
+namespace PHiLiP {
+namespace ODE {
+
+/// Runge-Kutta ODE solver (explicit or implicit) derived from ODESolver.
+#if PHILIP_DIM==1
+template <int dim, typename real, int n_rk_stages, typename MeshType = dealii::Triangulation<dim>>
+#else
+template <int dim, typename real, int n_rk_stages, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
+#endif
+class LowStorageRungeKuttaODESolver: public ODESolverBase <dim, real, MeshType>
+{
+public:
+    LowStorageRungeKuttaODESolver(std::shared_ptr< DGBase<dim, real, MeshType> > dg_input,
+            std::shared_ptr<LowStorageRKTableauBase<dim,real,MeshType>> rk_tableau_input,
+            std::shared_ptr<EmptyRRKBase<dim,real,MeshType>> RRK_object_input); ///< Constructor.
+
+    /// Function to evaluate solution update
+    double err_time_step(real dt, const bool pseudotime);
+
+    void step_in_time(real dt, const bool pseudotime);
+
+    /// Function to allocate the ODE system
+    void allocate_ode_system ();
+
+protected:
+    /// Stores Butcher tableau a and b, which specify the RK method
+    std::shared_ptr<LowStorageRKTableauBase<dim,real,MeshType>> butcher_tableau;
+
+    /// Stores functions related to relaxation Runge-Kutta (RRK).
+    /// Functions are empty by default.
+    std::shared_ptr<EmptyRRKBase<dim,real,MeshType>> relaxation_runge_kutta;
+
+    /// Implicit solver for diagonally-implicit RK methods, using Jacobian-free Newton-Krylov 
+    /** 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;
+
+    /// Storage for the derivative at each Runge-Kutta stage
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> rk_stage;
+
+    /// Indicator for zero diagonal elements; used to toggle implicit solve.
+    std::vector<bool> butcher_tableau_aii_is_zero;
+/*
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> storage_register_2;
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> storage_register_1;
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> storage_register_3;
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> storage_register_4;
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> rhs;
+*/ 
+
+    real w;
+
+    double epsilon[3];
+    int num_delta;
+    int rk_order;
+    double atol;
+    double rtol;
+    //double beta_controller[3];
+};
+
+} // ODE namespace
+} // PHiLiP namespace
+
+#endif
+

src/ode_solver/ode_solver_base.cpp

@@ -20,6 +20,15 @@ ODESolverBase<dim,real,MeshType>::ODESolverBase(std::shared_ptr< DGBase<dim, rea
         , pcout(std::cout, mpi_rank==0)
 {}
 
+
+template <int dim, typename real, typename MeshType> 
+double ODESolverBase<dim,real,MeshType>::err_time_step (real dt, const bool pseudotime)
+{
+    (void) pseudotime;
+    return dt;
+}
+
+
 template <int dim, typename real, typename MeshType>
 double ODESolverBase<dim,real,MeshType>::get_original_time_step() const
 {

src/ode_solver/ode_solver_base.h

@@ -42,6 +42,7 @@ class ODESolverBase
     /// Evaluate steady state solution.
     virtual int steady_state ();
 
+
     /// Ramps up the solution from p0 all the way up to the given global_final_poly_degree.
     /** This first interpolates the current solution to the P0 space as an initial solution.
      *  The P0 is then solved, interpolated to P1, and the process is repeated until the
@@ -63,6 +64,8 @@ class ODESolverBase
     /// Virtual function to evaluate solution update
     virtual void step_in_time(real dt, const bool pseudotime) = 0;
 
+    virtual double err_time_step (real dt, const bool pseudotime);
+
     /// Virtual function to allocate the ODE system
     virtual void allocate_ode_system () = 0;