From 4f19280a4b58dcc5ad9971e48d53de9bcc75e9be Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Wed, 12 Jun 2024 11:37:58 -0400 Subject: [PATCH 01/18] First commit of new Low Storage RK ODESolver. [Does not compile] --- .gitignore | 1 + src/ode_solver/CMakeLists.txt | 2 + .../low_storage_runge_kutta_ode_solver.cpp | 128 ++++++++++++++++++ .../low_storage_runge_kutta_ode_solver.h | 57 ++++++++ src/ode_solver/ode_solver_factory.cpp | 40 +++++- src/ode_solver/ode_solver_factory.h | 5 + .../low_storage_rk_tableau_base.cpp | 82 +++++++++++ .../low_storage_rk_tableau_base.h | 79 +++++++++++ src/parameters/parameters_ode_solver.cpp | 4 + src/parameters/parameters_ode_solver.h | 1 + 10 files changed, 395 insertions(+), 4 deletions(-) create mode 100644 src/ode_solver/low_storage_runge_kutta_ode_solver.cpp create mode 100644 src/ode_solver/low_storage_runge_kutta_ode_solver.h create mode 100644 src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp create mode 100644 src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h diff --git a/.gitignore b/.gitignore index 737a21955..4f25df95c 100644 --- a/.gitignore +++ b/.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 diff --git a/src/ode_solver/CMakeLists.txt b/src/ode_solver/CMakeLists.txt index 110b35231..5075df707 100644 --- a/src/ode_solver/CMakeLists.txt +++ b/src/ode_solver/CMakeLists.txt @@ -2,8 +2,10 @@ 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 relaxation_runge_kutta/empty_RRK_base.cpp relaxation_runge_kutta/runge_kutta_store_entropy.cpp relaxation_runge_kutta/rrk_ode_solver_base.cpp diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp new file mode 100644 index 000000000..1ed273d28 --- /dev/null +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -0,0 +1,128 @@ +#include "low_storage_runge_kutta_ode_solver.h" + +namespace PHiLiP { +namespace ODE { + +template +LowStorageRungeKuttaODESolver::LowStorageRungeKuttaODESolver(std::shared_ptr< DGBase > dg_input, + std::shared_ptr> rk_tableau_input, + std::shared_ptr> RRK_object_input) + : ODESolverBase(dg_input) + , butcher_tableau(rk_tableau_input) + , relaxation_runge_kutta(RRK_object_input) + , solver(dg_input) +{} + +template +void LowStorageRungeKuttaODESolver::step_in_time (real dt, const bool pseudotime) +{ + this->original_time_step = dt; + this->solution_update = this->dg->solution; //storing u_n + (void) pseudotime; + //My first change. + + //const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; + //double beta[6] = {0.0, 0.075152045700771, 0.211361016946069, 1.100713347634329, 0.728537814675568, 0.393172889823198}; + //double delta[7] = {1.0, 0.081252332929194, -1.083849060586449, -1.096110881845602, 2.859440022030827, -0.655568367959557, -0.194421504490852}; + // double beta_controller[3] = {0.70, -0.40, 0.0}; // PI34 + + dealii::LinearAlgebra::distributed::Vector s2; + s2.reinit(this->solution_update); + s2*=0; + dealii::LinearAlgebra::distributed::Vector s3; + s3.reinit(this->solution_update); + s3 = this->dg->solution; + dealii::LinearAlgebra::distributed::Vector s1 = s3; + //rhs.reinit(this->solution_update); + //rhs = s3; + dealii::LinearAlgebra::distributed::Vector u_hat = s2; + //u_hat.reinit(this->solution_update); + //u_hat = s2; + + dealii::LinearAlgebra::distributed::Vector rhs = s1; + + int m = 5; + double sum_delta = 0; + + for (int i = 1; i < m+1; i++ ){ + // s2 = s2 + delta[i-1] * s1; + s2.add(this->butcher_tableau->get_delta(i-1) , s1); + this->dg->solution = rhs; + this->dg->assemble_residual(); + this->dg->apply_inverse_global_mass_matrix(this->dg->right_hand_side, rhs); + // s1 = gamma[i][0] * s1 + gamma[i][1] * s2 + gamma[i][2] * s3 + beta[i] * dt * rhs; + s1 *= this->butcher_tableau->get_gamma(i, 0); + s1.add(this->butcher_tableau->get_gamma(i, 1), s2); + s1.add(this->butcher_tableau->get_gamma(i, 2), s3); + rhs *= dt; + s1.add(this->butcher_tableau->get_beta(i), rhs); + rhs = s1; + + } + for (int i = 0; ibutcher_tableau->get_delta(i); + } + // s2 = (s2 + delta[m] * s1 + delta[m+1] * s3) / sum_delta; + s2.add(this->butcher_tableau->get_delta(m), s1); + s2.add(this->butcher_tableau->get_delta(m+1), s3); + s2 /= sum_delta; + + this->dg->solution = s1; + u_hat = s2; + + ++(this->current_iteration); + this->current_time += dt; + this->pcout << "Time is: " << this->current_time < +void LowStorageRungeKuttaODESolver::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->rk_stage.resize(n_rk_stages); + for (int i=0; irk_stage[i].reinit(this->dg->solution); + } + + this->butcher_tableau->set_tableau(); + + this->butcher_tableau_aii_is_zero.resize(n_rk_stages); + std::fill(this->butcher_tableau_aii_is_zero.begin(), + this->butcher_tableau_aii_is_zero.end(), + false); + for (int i=0; ibutcher_tableau->get_a(i,i)==0.0) this->butcher_tableau_aii_is_zero[i] = true; + } +} + +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; +#if PHILIP_DIM != 1 + template class LowStorageRungeKuttaODESolver >; + template class LowStorageRungeKuttaODESolver >; + template class LowStorageRungeKuttaODESolver >; + template class LowStorageRungeKuttaODESolver >; +#endif + +} // ODESolver namespace +} // PHiLiP namespace diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.h b/src/ode_solver/low_storage_runge_kutta_ode_solver.h new file mode 100644 index 000000000..2a014d0ae --- /dev/null +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.h @@ -0,0 +1,57 @@ +#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 > +#else +template > +#endif +class LowStorageRungeKuttaODESolver: public ODESolverBase +{ +public: + LowStorageRungeKuttaODESolver(std::shared_ptr< DGBase > dg_input, + std::shared_ptr> rk_tableau_input, + std::shared_ptr> RRK_object_input); ///< Constructor. + + /// Function to evaluate solution update + 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> butcher_tableau; + + /// Stores functions related to relaxation Runge-Kutta (RRK). + /// Functions are empty by default. + std::shared_ptr> 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 solver; + + /// Storage for the derivative at each Runge-Kutta stage + std::vector> rk_stage; + + /// Indicator for zero diagonal elements; used to toggle implicit solve. + std::vector butcher_tableau_aii_is_zero; +}; + +} // ODE namespace +} // PHiLiP namespace + +#endif + diff --git a/src/ode_solver/ode_solver_factory.cpp b/src/ode_solver/ode_solver_factory.cpp index a012c960a..cb1ba9c5e 100644 --- a/src/ode_solver/ode_solver_factory.cpp +++ b/src/ode_solver/ode_solver_factory.cpp @@ -2,6 +2,7 @@ #include "parameters/all_parameters.h" #include "ode_solver_base.h" #include "runge_kutta_ode_solver.h" +#include "low_storage_runge_kutta_ode_solver.h" #include "implicit_ode_solver.h" #include "relaxation_runge_kutta/algebraic_rrk_ode_solver.h" #include "relaxation_runge_kutta/root_finding_rrk_ode_solver.h" @@ -10,6 +11,7 @@ #include #include "runge_kutta_methods/runge_kutta_methods.h" #include "runge_kutta_methods/rk_tableau_base.h" +#include "runge_kutta_methods/low_storage_rk_tableau_base.h" #include "relaxation_runge_kutta/empty_RRK_base.h" namespace PHiLiP { @@ -22,8 +24,8 @@ std::shared_ptr> ODESolverFactoryall_parameters->ode_solver_param.ode_solver_type; - if((ode_solver_type == ODEEnum::runge_kutta_solver)||(ode_solver_type == ODEEnum::rrk_explicit_solver)) - return create_RungeKuttaODESolver(dg_input); + if((ode_solver_type == ODEEnum::runge_kutta_solver)||(ode_solver_type == ODEEnum::rrk_explicit_solver || ode_solver_type == ODEEnum::low_storage_runge_kutta_solver)) + return create_RungeKuttaODESolver(dg_input); if(ode_solver_type == ODEEnum::implicit_solver) return std::make_shared>(dg_input); else { @@ -147,13 +149,43 @@ std::shared_ptr> ODESolverFactory> ls_rk_tableau = create_LowStorageRKTableau(dg_input);; + // ls_rk_tableau = std::make_shared> ("LSRK"); + // Hard-coded templating of n_rk_stages because it is not known at compile time + pcout << "Creating Low-Storage Runge Kutta ODE Solver with " + << n_rk_stages << " stage(s)..." << std::endl; + if (n_rk_stages == 1){ + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + } + else if (n_rk_stages == 2){ + + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + } + else if (n_rk_stages == 3){ + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + } + else if (n_rk_stages == 4){ + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + } + else{ + pcout << "Error: invalid number of stages. Aborting..." << std::endl; + std::abort(); + return nullptr; + } + } else { display_error_ode_solver_factory(ode_solver_type, false); return nullptr; } } +template +std::shared_ptr> ODESolverFactory::create_LowStorageRKTableau(std::shared_ptr< DGBase > dg_input) +{ + (void) dg_input; + return std::make_shared> ("LSRK"); +} + template std::shared_ptr> ODESolverFactory::create_RKTableau(std::shared_ptr< DGBase > dg_input) { diff --git a/src/ode_solver/ode_solver_factory.h b/src/ode_solver/ode_solver_factory.h index 90de00f25..517577c7b 100644 --- a/src/ode_solver/ode_solver_factory.h +++ b/src/ode_solver/ode_solver_factory.h @@ -6,6 +6,7 @@ #include "parameters/all_parameters.h" #include "reduced_order/pod_basis_base.h" #include "runge_kutta_methods/rk_tableau_base.h" +#include "runge_kutta_methods/low_storage_rk_tableau_base.h" #include "relaxation_runge_kutta/empty_RRK_base.h" namespace PHiLiP { @@ -42,6 +43,10 @@ class ODESolverFactory /// Creates an RKTableau object based on the specified RK method static std::shared_ptr> create_RKTableau(std::shared_ptr< DGBase > dg_input); + /// Creates an RKTableau object based on the specified RK method + static std::shared_ptr> create_LowStorageRKTableau(std::shared_ptr< DGBase > dg_input); + + /// Creates an RRK object with specified RRK type; if no RRK is being used, creates an RRK object with empty functions. /// Creates an RRK object with specified RRK type; if no RRK is being used, creates an RRK object with empty functions. static std::shared_ptr> create_RRKObject(std::shared_ptr< DGBase > dg_input, std::shared_ptr> rk_tableau); diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp new file mode 100644 index 000000000..d7db10bbf --- /dev/null +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp @@ -0,0 +1,82 @@ +#include "low_storage_rk_tableau_base.h" + + +namespace PHiLiP { +namespace ODE { + +template +LowStorageRKTableauBase :: LowStorageRKTableauBase ( + const std::string rk_method_string_input) + : rk_method_string(rk_method_string_input) + , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0) +{ + //this->butcher_tableau_a.reinit(n_rk_stages,n_rk_stages); + //this->butcher_tableau_b.reinit(n_rk_stages); + //this->butcher_tableau_c.reinit(n_rk_stages); + this->butcher_tableau_gamma.reinit(6,3); + this->butcher_tableau_beta.reinit(6); + this->butcher_tableau_delta.reinit(7); +} + +template +void LowStorageRKTableauBase :: set_tableau () +{ + set_gamma(); + set_beta(); + set_delta(); + pcout << "Assigned RK method: " << rk_method_string << std::endl; +} + +template +double LowStorageRKTableauBase :: get_gamma (const int i, const int j) const +{ + return butcher_tableau_gamma[i][j]; +} + +template +double LowStorageRKTableauBase :: get_beta (const int i) const +{ + return butcher_tableau_beta[i]; +} + +template +double LowStorageRKTableauBase :: get_delta (const int i) const +{ + return butcher_tableau_delta[i]; +} + +template +void LowStorageRKTableauBase :: set_gamma() +{ + const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; + for (int i = 0; i<6; i++){ + for (int j = 0; j<3; j++){ + this->butcher_tableau_gamma(i,j) = gamma[i][j]; + } + } +} + +template +void LowStorageRKTableauBase :: set_beta() +{ + double beta[6] = {0.0, 0.075152045700771, 0.211361016946069, 1.100713347634329, 0.728537814675568, 0.393172889823198}; + this->butcher_tableau_beta.fill(beta); +} + +template +void LowStorageRKTableauBase :: set_delta() +{ + double delta[7] = {1.0, 0.081252332929194, -1.083849060586449, -1.096110881845602, 2.859440022030827, -0.655568367959557, -0.194421504490852}; + this->butcher_tableau_delta.fill(delta); + +} + + +template class LowStorageRKTableauBase>; +template class LowStorageRKTableauBase>; +#if PHILIP_DIM != 1 +template class LowStorageRKTableauBase>; +#endif + +} // ODE namespace +} // PHiLiP namespace diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h new file mode 100644 index 000000000..b414340d0 --- /dev/null +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h @@ -0,0 +1,79 @@ +#ifndef __LOW_STORAGE_RK_TABLEAU_BASE__ +#define __LOW_STORAGE_RK_TABLEAU_BASE__ + +#include + +#include +#include +#include + +namespace PHiLiP { +namespace ODE { + +/// Base class for storing the RK method +#if PHILIP_DIM==1 +template > +#else +template > +#endif +class LowStorageRKTableauBase +{ +public: + /// Default constructor that will set the constants. + LowStorageRKTableauBase(const std::string rk_method_string_input); + + /// Destructor + virtual ~LowStorageRKTableauBase() = default; + +protected: + /// String identifying the RK method + const std::string rk_method_string; + + dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 + +public: + /// Returns Butcher tableau "a" coefficient at position [i][j] + double get_gamma(const int i, const int j) const; + + /// Returns Butcher tableau "b" coefficient at position [i] + double get_beta(const int i) const; + + /// Returns Butcher tableau "c" coefficient at position [i] + double get_delta(const int i) const; + + /// Calls setters for butcher tableau + void set_tableau(); + + +protected: + + + + /// Butcher tableau "a" + dealii::Table<2,double> butcher_tableau_gamma; + + /// Butcher tableau "b" + dealii::Table<1,double> butcher_tableau_beta; + + /// Butcher tableau "c" + dealii::Table<1,double> butcher_tableau_delta; + + /// Setter for butcher_tableau_a + // virtual void set_a() = 0; + void set_gamma(); + + /// Setter for butcher_tableau_b + // virtual void set_b() = 0; + void set_beta(); + + /// Setter for butcher_tableau_c + // virtual void set_c() = 0; + void set_delta(); + + +}; + +} // ODE namespace +} // PHiLiP namespace + +#endif diff --git a/src/parameters/parameters_ode_solver.cpp b/src/parameters/parameters_ode_solver.cpp index c7ca1cc87..5fda9056a 100644 --- a/src/parameters/parameters_ode_solver.cpp +++ b/src/parameters/parameters_ode_solver.cpp @@ -37,6 +37,7 @@ void ODESolverParam::declare_parameters (dealii::ParameterHandler &prm) prm.declare_entry("ode_solver_type", "implicit", dealii::Patterns::Selection( " runge_kutta | " + " low_storage_runge_kutta | " " implicit | " " rrk_explicit | " " pod_galerkin | " @@ -44,6 +45,7 @@ void ODESolverParam::declare_parameters (dealii::ParameterHandler &prm) "Type of ODE solver to use." "Choices are " " Date: Wed, 12 Jun 2024 12:28:45 -0400 Subject: [PATCH 02/18] low storage test --- ...me_refinement_study_advection_explicit.prm | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm diff --git a/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm b/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm new file mode 100644 index 000000000..41ecc15b3 --- /dev/null +++ b/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm @@ -0,0 +1,41 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions + +set dimension = 1 +set test_type = time_refinement_study +set pde_type = advection + +# Note: this was added to turn off check_same_coords() -- has no other function when dim!=1 +set use_periodic_bc = true + +# ODE solver +subsection ODE solver + set ode_solver_type = runge_kutta + set output_solution_every_dt_time_intervals = 0.1 + set initial_time_step = 5E-3 + set runge_kutta_method = low_storage_runge_kutta_solver +end + +subsection manufactured solution convergence study + # advection speed + set advection_0 = 1.0 + set advection_1 = 0.0 +end + + +subsection time_refinement_study + set number_of_times_to_solve = 5 + set refinement_ratio = 0.8 +end + +subsection flow_solver + set flow_case_type = periodic_1D_unsteady + set final_time = 1.0 + set poly_degree = 5 + subsection grid + set grid_left_bound = 0.0 + set grid_right_bound = 2.0 + set number_of_grid_elements_per_dimension = 32 + end +end From 01d72b34110101c726610de45c76b10e55c3c99f Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Wed, 12 Jun 2024 13:58:20 -0400 Subject: [PATCH 03/18] compiles --- .../low_storage_runge_kutta_ode_solver.cpp | 5 +++++ .../low_storage_rk_tableau_base.cpp | 4 ++-- .../time_refinement_study/CMakeLists.txt | 16 ++++++++++++++++ ..._time_refinement_study_advection_explicit.prm | 4 ++-- 4 files changed, 25 insertions(+), 4 deletions(-) diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 1ed273d28..fc4aa310a 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -91,8 +91,11 @@ void LowStorageRungeKuttaODESolver::allocate_ode_ this->dg->evaluate_mass_matrices(false); // creates and stores global mass matrix } } + this->pcout << std::endl; + this->butcher_tableau->set_tableau(); +/* this->rk_stage.resize(n_rk_stages); for (int i=0; irk_stage[i].reinit(this->dg->solution); @@ -106,7 +109,9 @@ void LowStorageRungeKuttaODESolver::allocate_ode_ false); for (int i=0; ibutcher_tableau->get_a(i,i)==0.0) this->butcher_tableau_aii_is_zero[i] = true; + } + */ } template class LowStorageRungeKuttaODESolver >; diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp index d7db10bbf..1253b4d3d 100644 --- a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp @@ -49,8 +49,8 @@ template void LowStorageRKTableauBase :: set_gamma() { const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; - for (int i = 0; i<6; i++){ - for (int j = 0; j<3; j++){ + for (int i = 0; i < 6; i++){ + for (int j = 0; j < 3; j++){ this->butcher_tableau_gamma(i,j) = gamma[i][j]; } } diff --git a/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt b/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt index 9f4f8f105..963311895 100644 --- a/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt +++ b/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt @@ -16,6 +16,22 @@ add_test( ) # ---------------------------------------- +# ======================================= +# Time Study (Linear Advection Low-Storage RK) +# ======================================= +# ---------------------------------------- +# Time refinement study on linear advection using a sinusoidal initial condition +# L2 error calculated with respect to the exact solution +# Test will fail if the convergence order is not close to the expected order +# ---------------------------------------- +configure_file(low_storage_rk_time_refinement_study_advection_explicit.prm low_storage_rk_time_refinement_study_advection_explicit.prm COPYONLY) +add_test( + NAME 1D_TIME_REFINEMENT_STUDY_ADVECTION_LOW_STORAGE + COMMAND mpirun -np 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/low_storage_rk_time_refinement_study_advection_explicit.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) +# ---------------------------------------- + # ======================================= # Time Study (Linear Advection Implicit RK) # ======================================= diff --git a/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm b/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm index 41ecc15b3..e735957bb 100644 --- a/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm +++ b/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm @@ -11,10 +11,10 @@ set use_periodic_bc = true # ODE solver subsection ODE solver - set ode_solver_type = runge_kutta + set ode_solver_type = low_storage_runge_kutta set output_solution_every_dt_time_intervals = 0.1 set initial_time_step = 5E-3 - set runge_kutta_method = low_storage_runge_kutta_solver + # set runge_kutta_method = low_storage_runge_kutta_solver end subsection manufactured solution convergence study From ca60d55c76d3528a7bf9c2d78a81b79d45c7e1d4 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Wed, 12 Jun 2024 14:17:04 -0400 Subject: [PATCH 04/18] this works --- .../runge_kutta_methods/low_storage_rk_tableau_base.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp index 1253b4d3d..4abd77579 100644 --- a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp @@ -51,7 +51,7 @@ void LowStorageRKTableauBase :: set_gamma() const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; for (int i = 0; i < 6; i++){ for (int j = 0; j < 3; j++){ - this->butcher_tableau_gamma(i,j) = gamma[i][j]; + this->butcher_tableau_gamma[i][j] = gamma[i][j]; } } } From 555565c1232808ac7882f94653bfe175920e04cc Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Thu, 13 Jun 2024 12:38:07 -0400 Subject: [PATCH 05/18] compiles --- .../low_storage_runge_kutta_ode_solver.cpp | 29 +++++++++++++++++-- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index fc4aa310a..16db49a19 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -19,7 +19,6 @@ void LowStorageRungeKuttaODESolver::step_in_time this->original_time_step = dt; this->solution_update = this->dg->solution; //storing u_n (void) pseudotime; - //My first change. //const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; //double beta[6] = {0.0, 0.075152045700771, 0.211361016946069, 1.100713347634329, 0.728537814675568, 0.393172889823198}; @@ -35,14 +34,23 @@ void LowStorageRungeKuttaODESolver::step_in_time dealii::LinearAlgebra::distributed::Vector s1 = s3; //rhs.reinit(this->solution_update); //rhs = s3; - dealii::LinearAlgebra::distributed::Vector u_hat = s2; + //dealii::LinearAlgebra::distributed::Vector u_hat = s2; //u_hat.reinit(this->solution_update); //u_hat = s2; dealii::LinearAlgebra::distributed::Vector rhs = s1; int m = 5; + int q_hat = 3; + int k = q_hat +1; double sum_delta = 0; + double error = 0.0; + double w = 0.0; + double atol = 0.001; + double rtol = 0.001; + double epsilon[3] = {1, 1, 1}; + double beta_controller[3] = {0.70, -0.40, 0}; + for (int i = 1; i < m+1; i++ ){ // s2 = s2 + delta[i-1] * s1; @@ -68,11 +76,26 @@ void LowStorageRungeKuttaODESolver::step_in_time s2 /= sum_delta; this->dg->solution = s1; - u_hat = s2; + //u_hat = s2; + + // error based step size + + for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { + error = s1.local_element(i) - s2.local_element(i); + w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); + } + w = pow(w / s1.local_size(), 1/2); + epsilon[2] = epsilon[1]; + epsilon[1] = epsilon[0]; + epsilon[0] = 1 / w; + dt = pow(epsilon[0], beta_controller[0]/k) * pow(epsilon[1], beta_controller[1]/k) * pow(epsilon[2], beta_controller[2]/k) * dt; + ++(this->current_iteration); this->current_time += dt; this->pcout << "Time is: " << this->current_time <pcout << "dt" << dt; + } From 8cb35d87f133bbd51ce85ae46beac6ad958647bc Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Thu, 13 Jun 2024 15:11:57 -0400 Subject: [PATCH 06/18] compiles but time step does not change as expected --- src/flow_solver/flow_solver.cpp | 6 +-- .../low_storage_runge_kutta_ode_solver.cpp | 44 +++++++++++++++---- .../low_storage_runge_kutta_ode_solver.h | 4 ++ src/ode_solver/ode_solver_base.cpp | 6 +++ src/ode_solver/ode_solver_base.h | 2 + src/ode_solver/runge_kutta_ode_solver.cpp | 3 +- 6 files changed, 53 insertions(+), 12 deletions(-) diff --git a/src/flow_solver/flow_solver.cpp b/src/flow_solver/flow_solver.cpp index 138b6ef57..15d616f6c 100644 --- a/src/flow_solver/flow_solver.cpp +++ b/src/flow_solver/flow_solver.cpp @@ -500,10 +500,10 @@ int FlowSolver::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); + next_time_step = ode_solver->err_time_step(time_step,false); // 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); diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 16db49a19..89c5c72f1 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -41,15 +41,17 @@ void LowStorageRungeKuttaODESolver::step_in_time dealii::LinearAlgebra::distributed::Vector rhs = s1; int m = 5; - int q_hat = 3; - int k = q_hat +1; + //int q_hat = 3; + //int k = q_hat +1; double sum_delta = 0; - double error = 0.0; - double w = 0.0; double atol = 0.001; double rtol = 0.001; - double epsilon[3] = {1, 1, 1}; - double beta_controller[3] = {0.70, -0.40, 0}; + double error = 0.0; + //double w = 0.0; + //double atol = 0.001; + //double rtol = 0.001; + //double epsilon[3] = {1, 1, 1}; + //double beta_controller[3] = {0.70, -0.40, 0}; for (int i = 1; i < m+1; i++ ){ @@ -80,6 +82,7 @@ void LowStorageRungeKuttaODESolver::step_in_time // error based step size + /* for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); @@ -88,8 +91,13 @@ void LowStorageRungeKuttaODESolver::step_in_time epsilon[2] = epsilon[1]; epsilon[1] = epsilon[0]; epsilon[0] = 1 / w; - dt = pow(epsilon[0], beta_controller[0]/k) * pow(epsilon[1], beta_controller[1]/k) * pow(epsilon[2], beta_controller[2]/k) * dt; - + */ + + for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { + error = s1.local_element(i) - s2.local_element(i); + w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); + } + w = pow(w / s1.local_size(), 1/2); ++(this->current_iteration); this->current_time += dt; @@ -97,7 +105,27 @@ void LowStorageRungeKuttaODESolver::step_in_time this->pcout << "dt" << dt; } +template +double LowStorageRungeKuttaODESolver::err_time_step (real dt, const bool pseudotime) +{ + (void) pseudotime; + int q_hat = 3; + int k = q_hat +1; + // double error = 0.0; + //double w = 0.0; + double epsilon[3] = {1, 1, 1}; + double beta_controller[3] = {0.70, -0.40, 0}; + // error based step size + + + epsilon[2] = epsilon[1]; + epsilon[1] = epsilon[0]; + epsilon[0] = 1 / w; + + dt = pow(epsilon[0], beta_controller[0]/k) * pow(epsilon[1], beta_controller[1]/k) * pow(epsilon[2], beta_controller[2]/k) * dt; + return dt; +} template void LowStorageRungeKuttaODESolver::allocate_ode_system () diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.h b/src/ode_solver/low_storage_runge_kutta_ode_solver.h index 2a014d0ae..80b37d1c9 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.h +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.h @@ -24,6 +24,8 @@ class LowStorageRungeKuttaODESolver: public ODESolverBase std::shared_ptr> 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 @@ -48,6 +50,8 @@ class LowStorageRungeKuttaODESolver: public ODESolverBase /// Indicator for zero diagonal elements; used to toggle implicit solve. std::vector butcher_tableau_aii_is_zero; + + real w; }; } // ODE namespace diff --git a/src/ode_solver/ode_solver_base.cpp b/src/ode_solver/ode_solver_base.cpp index 7353d521d..c4f6b403b 100644 --- a/src/ode_solver/ode_solver_base.cpp +++ b/src/ode_solver/ode_solver_base.cpp @@ -19,6 +19,12 @@ ODESolverBase::ODESolverBase(std::shared_ptr< DGBase +double ODESolverBase::err_time_step (real dt, const bool pseudotime) +{ + (void) pseudotime; + return dt; +} template double ODESolverBase::get_original_time_step() const diff --git a/src/ode_solver/ode_solver_base.h b/src/ode_solver/ode_solver_base.h index 92804b3e3..73f4eb992 100644 --- a/src/ode_solver/ode_solver_base.h +++ b/src/ode_solver/ode_solver_base.h @@ -42,6 +42,8 @@ class ODESolverBase /// Evaluate steady state solution. virtual int steady_state (); + virtual double err_time_step (real dt, const bool pseudotime); + /// 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 diff --git a/src/ode_solver/runge_kutta_ode_solver.cpp b/src/ode_solver/runge_kutta_ode_solver.cpp index 3943f5de9..6b0f65d90 100644 --- a/src/ode_solver/runge_kutta_ode_solver.cpp +++ b/src/ode_solver/runge_kutta_ode_solver.cpp @@ -151,7 +151,7 @@ void RungeKuttaODESolver::allocate_ode_system () } } this->pcout << std::endl; - + this->rk_stage.resize(n_rk_stages); for (int i=0; irk_stage[i].reinit(this->dg->solution); @@ -165,6 +165,7 @@ void RungeKuttaODESolver::allocate_ode_system () false); for (int i=0; ibutcher_tableau->get_a(i,i)==0.0) this->butcher_tableau_aii_is_zero[i] = true; + } } From 303cbf619e060b82cb37e2d61c32ddd431943147 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Thu, 13 Jun 2024 17:16:33 -0400 Subject: [PATCH 07/18] fixed w term need to fix dt --- .../low_storage_runge_kutta_ode_solver.cpp | 22 ++++++++++++++----- 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 89c5c72f1..dd45dfbc3 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -47,7 +47,7 @@ void LowStorageRungeKuttaODESolver::step_in_time double atol = 0.001; double rtol = 0.001; double error = 0.0; - //double w = 0.0; + w = 0.0; //double atol = 0.001; //double rtol = 0.001; //double epsilon[3] = {1, 1, 1}; @@ -95,14 +95,22 @@ void LowStorageRungeKuttaODESolver::step_in_time for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); + //this->pcout << "diff " << error; w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); + //this->pcout << "w " << w; } - w = pow(w / s1.local_size(), 1/2); - + this->pcout << std::endl; + //this->pcout << "w " << w; + w = pow(w / s1.local_size(), 0.5); + //this->pcout << " size s1 " << s1.local_size(); + //this->pcout << " math " << w/s1.local_size(); + //this->pcout << "w " << pow(w / s1.local_size(), 1/2); + //this->pcout << "w " << w; ++(this->current_iteration); this->current_time += dt; - this->pcout << "Time is: " << this->current_time <pcout << "dt" << dt; + this->pcout << " Time is: " << this->current_time <pcout << "w" << w; + //this->pcout << "dt" << dt; } template @@ -122,8 +130,10 @@ double LowStorageRungeKuttaODESolver::err_time_st epsilon[2] = epsilon[1]; epsilon[1] = epsilon[0]; epsilon[0] = 1 / w; - + //this->pcout << "dt" << dt; + //this->pcout << "epsilon" << epsilon[0]; dt = pow(epsilon[0], beta_controller[0]/k) * pow(epsilon[1], beta_controller[1]/k) * pow(epsilon[2], beta_controller[2]/k) * dt; + this->pcout << "dt" << dt; return dt; } From d858fe75b716d02c4d8375d3c49ef689513bf813 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Fri, 14 Jun 2024 16:15:33 -0400 Subject: [PATCH 08/18] this runs --- src/flow_solver/flow_solver.cpp | 6 +++++- .../low_storage_runge_kutta_ode_solver.cpp | 17 ++++++++++------- .../low_storage_runge_kutta_ode_solver.h | 2 ++ ...time_refinement_study_advection_explicit.prm | 4 ++-- 4 files changed, 19 insertions(+), 10 deletions(-) diff --git a/src/flow_solver/flow_solver.cpp b/src/flow_solver/flow_solver.cpp index 15d616f6c..829255bc6 100644 --- a/src/flow_solver/flow_solver.cpp +++ b/src/flow_solver/flow_solver.cpp @@ -501,18 +501,22 @@ int FlowSolver::run() const // update time step in flow_solver_case 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->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) { diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index dd45dfbc3..1187aea5d 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -11,7 +11,10 @@ LowStorageRungeKuttaODESolver::LowStorageRungeKu , 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 void LowStorageRungeKuttaODESolver::step_in_time (real dt, const bool pseudotime) @@ -44,8 +47,8 @@ void LowStorageRungeKuttaODESolver::step_in_time //int q_hat = 3; //int k = q_hat +1; double sum_delta = 0; - double atol = 0.001; - double rtol = 0.001; + double atol = 0.01; + double rtol = 0.01; double error = 0.0; w = 0.0; //double atol = 0.001; @@ -110,7 +113,7 @@ void LowStorageRungeKuttaODESolver::step_in_time this->current_time += dt; this->pcout << " Time is: " << this->current_time <pcout << "w" << w; - //this->pcout << "dt" << dt; + this->pcout << "dt" << dt; } template @@ -121,7 +124,6 @@ double LowStorageRungeKuttaODESolver::err_time_st int k = q_hat +1; // double error = 0.0; //double w = 0.0; - double epsilon[3] = {1, 1, 1}; double beta_controller[3] = {0.70, -0.40, 0}; // error based step size @@ -129,10 +131,11 @@ double LowStorageRungeKuttaODESolver::err_time_st epsilon[2] = epsilon[1]; epsilon[1] = epsilon[0]; - epsilon[0] = 1 / w; + epsilon[0] = 1.0 / w; //this->pcout << "dt" << dt; //this->pcout << "epsilon" << epsilon[0]; - dt = pow(epsilon[0], beta_controller[0]/k) * pow(epsilon[1], beta_controller[1]/k) * pow(epsilon[2], beta_controller[2]/k) * dt; + 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; + this->pcout << std::endl; this->pcout << "dt" << dt; return dt; } diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.h b/src/ode_solver/low_storage_runge_kutta_ode_solver.h index 80b37d1c9..83b210e6c 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.h +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.h @@ -52,6 +52,8 @@ class LowStorageRungeKuttaODESolver: public ODESolverBase std::vector butcher_tableau_aii_is_zero; real w; + + double epsilon[3]; }; } // ODE namespace diff --git a/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm b/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm index e735957bb..176b75277 100644 --- a/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm +++ b/tests/integration_tests_control_files/time_refinement_study/low_storage_rk_time_refinement_study_advection_explicit.prm @@ -13,7 +13,7 @@ set use_periodic_bc = true subsection ODE solver set ode_solver_type = low_storage_runge_kutta set output_solution_every_dt_time_intervals = 0.1 - set initial_time_step = 5E-3 + set initial_time_step = 1.0E-3 # set runge_kutta_method = low_storage_runge_kutta_solver end @@ -25,7 +25,7 @@ end subsection time_refinement_study - set number_of_times_to_solve = 5 + set number_of_times_to_solve = 1 set refinement_ratio = 0.8 end From 35170cc6235a25bf820ff601c9ddf8cbd2df6b0a Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Mon, 17 Jun 2024 11:51:04 -0400 Subject: [PATCH 09/18] step size change --- src/flow_solver/flow_solver.cpp | 2 +- .../low_storage_runge_kutta_ode_solver.cpp | 16 +++++++++++----- src/ode_solver/ode_solver_base.cpp | 3 +++ 3 files changed, 15 insertions(+), 6 deletions(-) diff --git a/src/flow_solver/flow_solver.cpp b/src/flow_solver/flow_solver.cpp index 829255bc6..a61d525da 100644 --- a/src/flow_solver/flow_solver.cpp +++ b/src/flow_solver/flow_solver.cpp @@ -502,7 +502,7 @@ int FlowSolver::run() const // update time step in flow_solver_case ode_solver->step_in_time(time_step,false); pcout << "time step " << time_step; - next_time_step = ode_solver->err_time_step(time_step,false); + next_time_step = ode_solver->err_time_step(time_step,false); // change back to err_time_step pcout << "next time step " << next_time_step; // advance solution // ode_solver->err_time_step(time_step); // pseudotime==false diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 1187aea5d..e2f7608e6 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -13,7 +13,7 @@ LowStorageRungeKuttaODESolver::LowStorageRungeKu , solver(dg_input) { epsilon[0] = 1.0; epsilon[1] = 1.0; - epsilon[2] = 1.0; + epsilon[2] = 1.0; } template @@ -46,11 +46,13 @@ void LowStorageRungeKuttaODESolver::step_in_time int m = 5; //int q_hat = 3; //int k = q_hat +1; + double sum_delta = 0; - double atol = 0.01; - double rtol = 0.01; + double atol = 0.1; + double rtol = 0.1; double error = 0.0; w = 0.0; + //double atol = 0.001; //double rtol = 0.001; //double epsilon[3] = {1, 1, 1}; @@ -85,7 +87,7 @@ void LowStorageRungeKuttaODESolver::step_in_time // error based step size - /* + for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); @@ -94,7 +96,7 @@ void LowStorageRungeKuttaODESolver::step_in_time epsilon[2] = epsilon[1]; epsilon[1] = epsilon[0]; epsilon[0] = 1 / w; - */ + for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); @@ -109,6 +111,7 @@ void LowStorageRungeKuttaODESolver::step_in_time //this->pcout << " math " << w/s1.local_size(); //this->pcout << "w " << pow(w / s1.local_size(), 1/2); //this->pcout << "w " << w; + ++(this->current_iteration); this->current_time += dt; this->pcout << " Time is: " << this->current_time <::step_in_time this->pcout << "dt" << dt; } + + template double LowStorageRungeKuttaODESolver::err_time_step (real dt, const bool pseudotime) { @@ -140,6 +145,7 @@ double LowStorageRungeKuttaODESolver::err_time_st return dt; } + template void LowStorageRungeKuttaODESolver::allocate_ode_system () { diff --git a/src/ode_solver/ode_solver_base.cpp b/src/ode_solver/ode_solver_base.cpp index c4f6b403b..03c8d59bb 100644 --- a/src/ode_solver/ode_solver_base.cpp +++ b/src/ode_solver/ode_solver_base.cpp @@ -19,6 +19,8 @@ ODESolverBase::ODESolverBase(std::shared_ptr< DGBase double ODESolverBase::err_time_step (real dt, const bool pseudotime) { @@ -26,6 +28,7 @@ double ODESolverBase::err_time_step (real dt, const bool pseu return dt; } + template double ODESolverBase::get_original_time_step() const { From 15ad87c04935e31dcc54860291aef666ed73e5ca Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Tue, 25 Jun 2024 10:33:42 -0400 Subject: [PATCH 10/18] save --- src/ode_solver/low_storage_runge_kutta_ode_solver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index e2f7608e6..feac9dc6a 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -87,7 +87,7 @@ void LowStorageRungeKuttaODESolver::step_in_time // error based step size - + /* for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); @@ -96,7 +96,7 @@ void LowStorageRungeKuttaODESolver::step_in_time epsilon[2] = epsilon[1]; epsilon[1] = epsilon[0]; epsilon[0] = 1 / w; - + */ for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); From 820d80c0e8f092c9bfd92bcb20968b8a51dd3f9f Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Tue, 25 Jun 2024 14:11:47 -0400 Subject: [PATCH 11/18] works with mpi > 1 --- .../low_storage_runge_kutta_ode_solver.cpp | 16 +++++++++++----- ...D_inviscid_isentropic_vortex_h_refinement.prm | 4 ++-- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index feac9dc6a..f083761d1 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -48,8 +48,8 @@ void LowStorageRungeKuttaODESolver::step_in_time //int k = q_hat +1; double sum_delta = 0; - double atol = 0.1; - double rtol = 0.1; + double atol = 0.001; + double rtol = 0.001; double error = 0.0; w = 0.0; @@ -97,20 +97,26 @@ void LowStorageRungeKuttaODESolver::step_in_time epsilon[1] = epsilon[0]; epsilon[0] = 1 / w; */ - + double global_size = dealii::Utilities::MPI::sum(s1.local_size(), this->mpi_communicator); for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { error = s1.local_element(i) - s2.local_element(i); //this->pcout << "diff " << error; w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); //this->pcout << "w " << w; } + // loop sums elements at each mpi processor this->pcout << std::endl; //this->pcout << "w " << w; - w = pow(w / s1.local_size(), 0.5); + + w = dealii::Utilities::MPI::sum(w, this->mpi_communicator); + // sums over all elements + w = pow(w / global_size, 0.5); + //this->pcout << " size s1 " << s1.local_size(); //this->pcout << " math " << w/s1.local_size(); //this->pcout << "w " << pow(w / s1.local_size(), 1/2); - //this->pcout << "w " << w; + this->pcout << std::endl; + std::cout << "w " << w; ++(this->current_iteration); this->current_time += dt; diff --git a/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm b/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm index b4d7f3f94..129130330 100644 --- a/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm +++ b/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm @@ -25,8 +25,8 @@ set two_point_num_flux_type = IR # ODE solver subsection ODE solver set ode_output = quiet - set ode_solver_type = runge_kutta - set runge_kutta_method = rk4_ex + set ode_solver_type = low_storage_runge_kutta +# set runge_kutta_method = rk4_ex end subsection linear solver From e758b4a2e8012036577857e0148621e3dc48ca11 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Wed, 26 Jun 2024 13:44:06 -0400 Subject: [PATCH 12/18] update --- src/flow_solver/flow_solver.cpp | 2 +- .../low_storage_runge_kutta_ode_solver.cpp | 127 +++++------------- src/ode_solver/ode_solver_base.h | 3 +- .../3D_navier_stokes_RRK_check.prm | 6 +- 4 files changed, 41 insertions(+), 97 deletions(-) diff --git a/src/flow_solver/flow_solver.cpp b/src/flow_solver/flow_solver.cpp index a61d525da..3ee9fb67b 100644 --- a/src/flow_solver/flow_solver.cpp +++ b/src/flow_solver/flow_solver.cpp @@ -502,7 +502,7 @@ int FlowSolver::run() const // update time step in flow_solver_case ode_solver->step_in_time(time_step,false); pcout << "time step " << time_step; - next_time_step = ode_solver->err_time_step(time_step,false); // change back to err_time_step + next_time_step = ode_solver->err_time_step(time_step,false); pcout << "next time step " << next_time_step; // advance solution // ode_solver->err_time_step(time_step); // pseudotime==false diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index f083761d1..32834fb61 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -23,107 +23,71 @@ void LowStorageRungeKuttaODESolver::step_in_time this->solution_update = this->dg->solution; //storing u_n (void) pseudotime; - //const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; - //double beta[6] = {0.0, 0.075152045700771, 0.211361016946069, 1.100713347634329, 0.728537814675568, 0.393172889823198}; - //double delta[7] = {1.0, 0.081252332929194, -1.083849060586449, -1.096110881845602, 2.859440022030827, -0.655568367959557, -0.194421504490852}; - // double beta_controller[3] = {0.70, -0.40, 0.0}; // PI34 - - dealii::LinearAlgebra::distributed::Vector s2; - s2.reinit(this->solution_update); - s2*=0; - dealii::LinearAlgebra::distributed::Vector s3; - s3.reinit(this->solution_update); - s3 = this->dg->solution; - dealii::LinearAlgebra::distributed::Vector s1 = s3; - //rhs.reinit(this->solution_update); - //rhs = s3; - //dealii::LinearAlgebra::distributed::Vector u_hat = s2; - //u_hat.reinit(this->solution_update); - //u_hat = s2; - - dealii::LinearAlgebra::distributed::Vector rhs = s1; - + dealii::LinearAlgebra::distributed::Vector storage_register_2; + storage_register_2.reinit(this->solution_update); + storage_register_2*=0; + dealii::LinearAlgebra::distributed::Vector storage_register_3; + storage_register_3.reinit(this->solution_update); + storage_register_3 = this->dg->solution; + dealii::LinearAlgebra::distributed::Vector storage_register_1 = storage_register_3; + + dealii::LinearAlgebra::distributed::Vector rhs = storage_register_1; + int m = 5; - //int q_hat = 3; - //int k = q_hat +1; double sum_delta = 0; double atol = 0.001; double rtol = 0.001; double error = 0.0; w = 0.0; - - //double atol = 0.001; - //double rtol = 0.001; - //double epsilon[3] = {1, 1, 1}; - //double beta_controller[3] = {0.70, -0.40, 0}; - for (int i = 1; i < m+1; i++ ){ - // s2 = s2 + delta[i-1] * s1; - s2.add(this->butcher_tableau->get_delta(i-1) , s1); + // 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); - // s1 = gamma[i][0] * s1 + gamma[i][1] * s2 + gamma[i][2] * s3 + beta[i] * dt * rhs; - s1 *= this->butcher_tableau->get_gamma(i, 0); - s1.add(this->butcher_tableau->get_gamma(i, 1), s2); - s1.add(this->butcher_tableau->get_gamma(i, 2), s3); + // 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; - s1.add(this->butcher_tableau->get_beta(i), rhs); - rhs = s1; + storage_register_1.add(this->butcher_tableau->get_beta(i), rhs); + rhs = storage_register_1; } - for (int i = 0; ibutcher_tableau->get_delta(i); } - // s2 = (s2 + delta[m] * s1 + delta[m+1] * s3) / sum_delta; - s2.add(this->butcher_tableau->get_delta(m), s1); - s2.add(this->butcher_tableau->get_delta(m+1), s3); - s2 /= sum_delta; - - this->dg->solution = s1; - //u_hat = s2; + // storage_register_2 = (storage_register_2 + delta[m] * storage_register_1 + delta[m+1] * storage_register_3) / sum_delta; + storage_register_2.add(this->butcher_tableau->get_delta(m), storage_register_1); + storage_register_2.add(this->butcher_tableau->get_delta(m+1), storage_register_3); + storage_register_2 /= sum_delta; - // error based step size + this->dg->solution = storage_register_1; + + double global_size = dealii::Utilities::MPI::sum(storage_register_1.local_size(), this->mpi_communicator); - /* - for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { - error = s1.local_element(i) - s2.local_element(i); - w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); - } - w = pow(w / s1.local_size(), 1/2); - epsilon[2] = epsilon[1]; - epsilon[1] = epsilon[0]; - epsilon[0] = 1 / w; - */ - double global_size = dealii::Utilities::MPI::sum(s1.local_size(), this->mpi_communicator); - for (dealii::LinearAlgebra::distributed::Vector::size_type i = 0; i < s1.local_size(); ++i) { - error = s1.local_element(i) - s2.local_element(i); - //this->pcout << "diff " << error; - w = w + pow(error / (atol + rtol * std::max(std::abs(s1.local_element(i)), std::abs(s2.local_element(i)))), 2); - //this->pcout << "w " << w; - } // loop sums elements at each mpi processor + for (dealii::LinearAlgebra::distributed::Vector::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); + } + this->pcout << std::endl; - //this->pcout << "w " << w; + // sum over all elements w = dealii::Utilities::MPI::sum(w, this->mpi_communicator); - // sums over all elements + w = pow(w / global_size, 0.5); - //this->pcout << " size s1 " << s1.local_size(); - //this->pcout << " math " << w/s1.local_size(); - //this->pcout << "w " << pow(w / s1.local_size(), 1/2); this->pcout << std::endl; - std::cout << "w " << w; + // std::cout << "w " << w; ++(this->current_iteration); this->current_time += dt; this->pcout << " Time is: " << this->current_time <pcout << "w" << w; this->pcout << "dt" << dt; - } @@ -133,8 +97,6 @@ double LowStorageRungeKuttaODESolver::err_time_st (void) pseudotime; int q_hat = 3; int k = q_hat +1; - // double error = 0.0; - //double w = 0.0; double beta_controller[3] = {0.70, -0.40, 0}; // error based step size @@ -143,11 +105,9 @@ double LowStorageRungeKuttaODESolver::err_time_st epsilon[2] = epsilon[1]; epsilon[1] = epsilon[0]; epsilon[0] = 1.0 / w; - //this->pcout << "dt" << dt; - //this->pcout << "epsilon" << epsilon[0]; 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; this->pcout << std::endl; - this->pcout << "dt" << dt; + //this->pcout << "dt" << dt; return dt; } @@ -171,23 +131,6 @@ void LowStorageRungeKuttaODESolver::allocate_ode_ this->pcout << std::endl; this->butcher_tableau->set_tableau(); -/* - this->rk_stage.resize(n_rk_stages); - for (int i=0; irk_stage[i].reinit(this->dg->solution); - } - - this->butcher_tableau->set_tableau(); - - this->butcher_tableau_aii_is_zero.resize(n_rk_stages); - std::fill(this->butcher_tableau_aii_is_zero.begin(), - this->butcher_tableau_aii_is_zero.end(), - false); - for (int i=0; ibutcher_tableau->get_a(i,i)==0.0) this->butcher_tableau_aii_is_zero[i] = true; - - } - */ } template class LowStorageRungeKuttaODESolver >; diff --git a/src/ode_solver/ode_solver_base.h b/src/ode_solver/ode_solver_base.h index 73f4eb992..bea920ee6 100644 --- a/src/ode_solver/ode_solver_base.h +++ b/src/ode_solver/ode_solver_base.h @@ -42,7 +42,6 @@ class ODESolverBase /// Evaluate steady state solution. virtual int steady_state (); - virtual double err_time_step (real dt, const bool pseudotime); /// 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. @@ -65,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; diff --git a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/3D_navier_stokes_RRK_check.prm b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/3D_navier_stokes_RRK_check.prm index c3f7327ef..4d1cfff8f 100644 --- a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/3D_navier_stokes_RRK_check.prm +++ b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/3D_navier_stokes_RRK_check.prm @@ -29,8 +29,8 @@ set diss_num_flux = symm_internal_penalty # ODE solver subsection ODE solver set ode_output = quiet - set ode_solver_type = rrk_explicit - set runge_kutta_method = ssprk3_ex + set ode_solver_type = low_storage_runge_kutta + #set runge_kutta_method = ssprk3_ex subsection rrk root solver set rrk_root_solver_output = verbose end @@ -54,7 +54,7 @@ end subsection flow_solver set flow_case_type = taylor_green_vortex set poly_degree = 2 - set final_time = 0.1 + set final_time = 10 set courant_friedrichs_lewy_number = 0.1 set unsteady_data_table_filename = turbulent_quantities subsection grid From ee073a4c8b1c31e85dd69b41a9bca596db08bce0 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Thu, 27 Jun 2024 13:37:09 -0400 Subject: [PATCH 13/18] add test euler inviscid tgv largest stable cfl --- .../CMakeLists.txt | 14 +++++ .../euler_inviscid_TGV_largest_stable_CFL.prm | 57 +++++++++++++++++++ 2 files changed, 71 insertions(+) create mode 100644 tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm diff --git a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/CMakeLists.txt b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/CMakeLists.txt index dd00dfdc3..af79be78d 100644 --- a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/CMakeLists.txt +++ b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/CMakeLists.txt @@ -1,3 +1,17 @@ +# ======================================= +# Euler Inviscid TGV Largest Stable CFL +# ======================================= +# ---------------------------------------- +# +# ---------------------------------------- + +configure_file(euler_inviscid_TGV_largest_stable_CFL.prm euler_inviscid_TGV_largest_stable_CFL.prm COPYONLY) +add_test( + NAME MPI_3D_EULER_INVISCID_TGV_LARGEST_STABLE_CFL + COMMAND mpirun -n ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/euler_inviscid_TGV_largest_stable_CFL.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + # ======================================= # Energy conservation test for inviscid Burgers # ======================================= diff --git a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm new file mode 100644 index 000000000..dc64e5491 --- /dev/null +++ b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm @@ -0,0 +1,57 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions + +set dimension = 3 +set run_type = flow_simulation +set pde_type = euler + +# DG formulation +set use_weak_form = false +set use_split_form = true +set flux_nodes_type = GLL + +set flux_reconstruction = cDG +set use_inverse_mass_on_the_fly = true + +# Note: this was added to turn off check_same_coords() -- has no other function when dim!=1 +set use_periodic_bc = true + +# numerical fluxes +set conv_num_flux = two_point_flux +set two_point_num_flux_type = Ra + +# ODE solver +subsection ODE solver + set ode_output = quiet + set ode_solver_type = low_storage_runge_kutta + #set ode_solver_type = rrk_explicit + #set runge_kutta_method = ssprk3_ex + subsection rrk root solver + set rrk_root_solver_output = verbose + end +end + +# freestream Mach number +subsection euler + set mach_infinity = 0.1 +end + +subsection flow_solver + set flow_case_type = taylor_green_vortex + set poly_degree = 3 + set final_time = 14 + set courant_friedrichs_lewy_number = 0.48 + set unsteady_data_table_filename = tgv_unsteady + set adaptive_time_step = true + set end_exactly_at_final_time = false + subsection grid + set grid_left_bound = 0.0 + set grid_right_bound = 6.28318530717958623200 + set number_of_grid_elements_per_dimension = 8 + end + subsection taylor_green_vortex + set do_calculate_numerical_entropy = true + set density_initial_condition_type = isothermal + end +end From 7acd2543d1e585db389c71d9ec41731f0c9d9029 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Tue, 2 Jul 2024 14:03:28 -0400 Subject: [PATCH 14/18] add 3s*+ method and RK3(2)5F[3S*+] --- src/ode_solver/CMakeLists.txt | 1 + .../low_storage_runge_kutta_ode_solver.cpp | 57 +++++--- src/ode_solver/ode_solver_factory.cpp | 18 ++- .../low_storage_rk_tableau_base.cpp | 20 ++- .../low_storage_rk_tableau_base.h | 15 ++- .../low_storage_runge_kutta_methods.cpp | 126 ++++++++++++++++++ .../low_storage_runge_kutta_methods.h | 54 ++++++++ src/parameters/parameters_ode_solver.cpp | 11 ++ src/parameters/parameters_ode_solver.h | 5 +- 9 files changed, 280 insertions(+), 27 deletions(-) create mode 100644 src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp create mode 100644 src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.h diff --git a/src/ode_solver/CMakeLists.txt b/src/ode_solver/CMakeLists.txt index 5075df707..5f864d51a 100644 --- a/src/ode_solver/CMakeLists.txt +++ b/src/ode_solver/CMakeLists.txt @@ -6,6 +6,7 @@ set(ODE_SOURCE 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 diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 32834fb61..38b5dd1de 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -32,8 +32,16 @@ void LowStorageRungeKuttaODESolver::step_in_time dealii::LinearAlgebra::distributed::Vector storage_register_1 = storage_register_3; dealii::LinearAlgebra::distributed::Vector rhs = storage_register_1; + + dealii::LinearAlgebra::distributed::Vector storage_register_4; + + if (this->butcher_tableau->get_b_hat(n_rk_stages+1) != 0){ + storage_register_4.reinit(this->solution_update); + storage_register_4 = this->dg->solution; + } + - int m = 5; + // int m = 5; double sum_delta = 0; double atol = 0.001; @@ -41,7 +49,7 @@ void LowStorageRungeKuttaODESolver::step_in_time double error = 0.0; w = 0.0; - for (int i = 1; i < m+1; i++ ){ + 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; @@ -53,27 +61,46 @@ void LowStorageRungeKuttaODESolver::step_in_time 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(n_rk_stages + 1) != 0){ + storage_register_4.add(this->butcher_tableau->get_b_hat(i), rhs); + } rhs = storage_register_1; } - for (int i = 0; i < m+2; i++){ + + // 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 < n_rk_stages +2; i++){ 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 + } + else if (this->butcher_tableau->get_b_hat(n_rk_stages+1) != 0){ + storage_register_1 *= dt; + storage_register_4.add(this->butcher_tableau->get_b_hat(n_rk_stages + 1), storage_register_1); + // u_hat = s1 } - // storage_register_2 = (storage_register_2 + delta[m] * storage_register_1 + delta[m+1] * storage_register_3) / sum_delta; - storage_register_2.add(this->butcher_tableau->get_delta(m), storage_register_1); - storage_register_2.add(this->butcher_tableau->get_delta(m+1), storage_register_3); - storage_register_2 /= sum_delta; - this->dg->solution = storage_register_1; double global_size = dealii::Utilities::MPI::sum(storage_register_1.local_size(), this->mpi_communicator); - // loop sums elements at each mpi processor - for (dealii::LinearAlgebra::distributed::Vector::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); + if (this->butcher_tableau->get_b_hat(1) == 0){ + // loop sums elements at each mpi processor + for (dealii::LinearAlgebra::distributed::Vector::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::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 @@ -92,12 +119,12 @@ void LowStorageRungeKuttaODESolver::step_in_time template -double LowStorageRungeKuttaODESolver::err_time_step (real dt, const bool pseudotime) +double LowStorageRungeKuttaODESolver::err_time_step (real dt, const bool pseudotime) { (void) pseudotime; int q_hat = 3; int k = q_hat +1; - double beta_controller[3] = {0.70, -0.40, 0}; + double beta_controller[3] = {0.70, -0.23, 0}; // error based step size diff --git a/src/ode_solver/ode_solver_factory.cpp b/src/ode_solver/ode_solver_factory.cpp index cb1ba9c5e..f2d6b77bf 100644 --- a/src/ode_solver/ode_solver_factory.cpp +++ b/src/ode_solver/ode_solver_factory.cpp @@ -12,6 +12,7 @@ #include "runge_kutta_methods/runge_kutta_methods.h" #include "runge_kutta_methods/rk_tableau_base.h" #include "runge_kutta_methods/low_storage_rk_tableau_base.h" +#include "runge_kutta_methods/low_storage_runge_kutta_methods.h" #include "relaxation_runge_kutta/empty_RRK_base.h" namespace PHiLiP { @@ -182,8 +183,21 @@ std::shared_ptr> ODESolverFactory std::shared_ptr> ODESolverFactory::create_LowStorageRKTableau(std::shared_ptr< DGBase > dg_input) { - (void) dg_input; - return std::make_shared> ("LSRK"); + // (void) dg_input; + dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0); + using RKMethodEnum = Parameters::ODESolverParam::RKMethodEnum; + const RKMethodEnum rk_method = dg_input->all_parameters->ode_solver_param.runge_kutta_method; + + const int n_rk_stages = dg_input->all_parameters->ode_solver_param.n_rk_stages; + const int num_delta = dg_input->all_parameters->ode_solver_param.num_delta; + if (rk_method == RKMethodEnum::RK3_2_5F_3SStarPlus) return std::make_shared> (n_rk_stages, num_delta, "RK3_2_5F_3SStarPlus"); + if (rk_method == RKMethodEnum::RK4_3_5_3SStar) return std::make_shared> (n_rk_stages, num_delta, "RK4_3_5_3SStar"); + //return std::make_shared> ("LSRK"); + else { + pcout << "Error: invalid RK method. Aborting..." << std::endl; + std::abort(); + return nullptr; + } } template diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp index 4abd77579..b358e04a0 100644 --- a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp @@ -5,7 +5,7 @@ namespace PHiLiP { namespace ODE { template -LowStorageRKTableauBase :: LowStorageRKTableauBase ( +LowStorageRKTableauBase :: LowStorageRKTableauBase (const int n_rk_stages, const int num_delta, const std::string rk_method_string_input) : rk_method_string(rk_method_string_input) , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0) @@ -13,9 +13,10 @@ LowStorageRKTableauBase :: LowStorageRKTableauBase ( //this->butcher_tableau_a.reinit(n_rk_stages,n_rk_stages); //this->butcher_tableau_b.reinit(n_rk_stages); //this->butcher_tableau_c.reinit(n_rk_stages); - this->butcher_tableau_gamma.reinit(6,3); - this->butcher_tableau_beta.reinit(6); - this->butcher_tableau_delta.reinit(7); + this->butcher_tableau_gamma.reinit(n_rk_stages+1,3); + this->butcher_tableau_beta.reinit(n_rk_stages+1); + this->butcher_tableau_delta.reinit(num_delta); + this->butcher_tableau_b_hat.reinit(n_rk_stages+2); } template @@ -45,9 +46,18 @@ double LowStorageRKTableauBase :: get_delta (const int i) co return butcher_tableau_delta[i]; } +template +double LowStorageRKTableauBase :: get_b_hat (const int i) const +{ + return butcher_tableau_b_hat[i]; +} + + +/* template void LowStorageRKTableauBase :: set_gamma() { + if const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; for (int i = 0; i < 6; i++){ for (int j = 0; j < 3; j++){ @@ -70,7 +80,7 @@ void LowStorageRKTableauBase :: set_delta() this->butcher_tableau_delta.fill(delta); } - +*/ template class LowStorageRKTableauBase>; template class LowStorageRKTableauBase>; diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h index b414340d0..5ebc89c1b 100644 --- a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.h @@ -20,7 +20,7 @@ class LowStorageRKTableauBase { public: /// Default constructor that will set the constants. - LowStorageRKTableauBase(const std::string rk_method_string_input); + LowStorageRKTableauBase(const int n_rk_stages, const int num_delta, const std::string rk_method_string_input); /// Destructor virtual ~LowStorageRKTableauBase() = default; @@ -41,6 +41,8 @@ class LowStorageRKTableauBase /// Returns Butcher tableau "c" coefficient at position [i] double get_delta(const int i) const; + double get_b_hat(const int i) const; + /// Calls setters for butcher tableau void set_tableau(); @@ -57,18 +59,23 @@ class LowStorageRKTableauBase /// Butcher tableau "c" dealii::Table<1,double> butcher_tableau_delta; + + /// Butcher tableau "c" + dealii::Table<1,double> butcher_tableau_b_hat; /// Setter for butcher_tableau_a // virtual void set_a() = 0; - void set_gamma(); + virtual void set_gamma() = 0; /// Setter for butcher_tableau_b // virtual void set_b() = 0; - void set_beta(); + virtual void set_beta() = 0; /// Setter for butcher_tableau_c // virtual void set_c() = 0; - void set_delta(); + virtual void set_delta() = 0; + + virtual void set_b_hat() = 0; }; diff --git a/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp b/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp new file mode 100644 index 000000000..56424b1b6 --- /dev/null +++ b/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp @@ -0,0 +1,126 @@ +#include "low_storage_runge_kutta_methods.h" + +namespace PHiLiP { +namespace ODE { + +//################################################################## + +// RK4(3)5[3S*] + +template +void RK4_3_5_3SStar :: set_gamma() +{ + const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; + for (int i = 0; i < 6; i++){ + for (int j = 0; j < 3; j++){ + this->butcher_tableau_gamma[i][j] = gamma[i][j]; + } + } +} + +template +void RK4_3_5_3SStar :: set_beta() +{ + const double beta[6] = {0.0, 0.075152045700771, 0.211361016946069, 1.100713347634329, 0.728537814675568, 0.393172889823198}; + this->butcher_tableau_beta.fill(beta); +} + +template +void RK4_3_5_3SStar :: set_delta() +{ + const double delta[7] = {1.0, 0.081252332929194, -1.083849060586449, -1.096110881845602, 2.859440022030827, -0.655568367959557, -0.194421504490852}; + this->butcher_tableau_delta.fill(delta); + +} + +template +void RK4_3_5_3SStar :: set_b_hat() +{ + const double b_hat[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + this->butcher_tableau_b_hat.fill(b_hat); +} + +//################################################################## + +// RK3(2)5F[3S*+] + +template +void RK3_2_5F_3SStarPlus :: set_gamma() +{ + const double gamma[5][3] = {{0.0, 1.0, 0.0}, {0.2587771979725733308135192812685323706, 0.5528354909301389892439698870483746541, 0.0}, {-0.1324380360140723382965420909764953437, 0.6731871608203061824849561782794643600, 0.0}, {0.05056033948190826045833606441415585735, 0.2803103963297672407841316576323901761, 0.2752563273304676380891217287572780582}, {0.5670532000739313812633197158607642990, 0.5521525447020610386070346724931300367, -0.8950526174674033822276061734289327568}}; + for (int i = 0; i < 5; i++){ + for (int j = 0; j < 3; j++){ + this->butcher_tableau_gamma[i][j] = gamma[i][j]; + } + } +} + +template +void RK3_2_5F_3SStarPlus :: set_beta() +{ + const double beta[5] = {0.11479359710235412, 0.089334428531133159, 0.43558710250086169, 0.24735761882014512, 0.058819144221557401}; + this->butcher_tableau_beta.fill(beta); +} + +template +void RK3_2_5F_3SStarPlus :: set_delta() +{ + const double delta[5] = {1.0, 0.34076558793345252, 0.34143826550033862, 0.72292753667879872, 0.0}; + this->butcher_tableau_delta.fill(delta); + +} + +template +void RK3_2_5F_3SStarPlus :: set_b_hat() +{ + const double b_hat[6] = {0.094841667050357029, 0.17263713394303537, 0.39982431890843712, 0.17180168075801786, 0.058819144221557401, 0.1020760551185952388626787099944507877}; + this->butcher_tableau_b_hat.fill(b_hat); +} + +//################################################################## +template class RK4_3_5_3SStar >; +template class RK4_3_5_3SStar >; +#if PHILIP_DIM != 1 + template class RK4_3_5_3SStar >; +#endif + +template class RK3_2_5F_3SStarPlus >; +template class RK3_2_5F_3SStarPlus >; +#if PHILIP_DIM != 1 + template class RK3_2_5F_3SStarPlus >; +#endif + +/* +template class EulerExplicit >; +template class EulerExplicit >; +#if PHILIP_DIM != 1 + template class EulerExplicit >; +#endif + +template class HeunExplicit >; +template class HeunExplicit >; +#if PHILIP_DIM != 1 + template class HeunExplicit >; +#endif + +template class EulerImplicit >; +template class EulerImplicit >; +#if PHILIP_DIM != 1 + template class EulerImplicit >; +#endif + +template class DIRK2Implicit >; +template class DIRK2Implicit >; +#if PHILIP_DIM != 1 + template class DIRK2Implicit >; +#endif + +template class DIRK3Implicit >; +template class DIRK3Implicit >; +#if PHILIP_DIM != 1 + template class DIRK3Implicit >; +#endif +*/ + +} // ODESolver namespace +} // PHiLiP namespace diff --git a/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.h b/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.h new file mode 100644 index 000000000..6bb1b9aa4 --- /dev/null +++ b/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.h @@ -0,0 +1,54 @@ +#ifndef __LOW_STORAGE_RUNGE_KUTTA_METHODS__ +#define __LOW_STORAGE_RUNGE_KUTTA_METHODS__ + +#include "low_storage_rk_tableau_base.h" + +namespace PHiLiP { +namespace ODE { + +#if PHILIP_DIM==1 +template > +#else +template > +#endif +class RK4_3_5_3SStar: public LowStorageRKTableauBase +{ +public: + /// Constructor + RK4_3_5_3SStar(const int n_rk_stages, const int num_delta, const std::string rk_method_string_input) + : LowStorageRKTableauBase(n_rk_stages, num_delta, rk_method_string_input) { } + +protected: + /// Setters + void set_gamma() override; + void set_beta() override; + void set_delta() override; + void set_b_hat() override; +}; + + +#if PHILIP_DIM==1 +template > +#else +template > +#endif +class RK3_2_5F_3SStarPlus: public LowStorageRKTableauBase +{ +public: + /// Constructor + RK3_2_5F_3SStarPlus(const int n_rk_stages, const int num_delta, const std::string rk_method_string_input) + : LowStorageRKTableauBase(n_rk_stages, num_delta, rk_method_string_input) { } + +protected: + /// Setters + void set_gamma() override; + void set_beta() override; + void set_delta() override; + void set_b_hat() override; +}; + + +} // ODE namespace +} // PHiLiP namespace + +#endif diff --git a/src/parameters/parameters_ode_solver.cpp b/src/parameters/parameters_ode_solver.cpp index 5fda9056a..76b2fd23f 100644 --- a/src/parameters/parameters_ode_solver.cpp +++ b/src/parameters/parameters_ode_solver.cpp @@ -213,6 +213,17 @@ void ODESolverParam::parse_parameters (dealii::ParameterHandler &prm) n_rk_stages = 3; rk_order = 3; } + else if (rk_method_string == "RK3_2_5F_3SStarPlus"){ + runge_kutta_method = RKMethodEnum::RK3_2_5F_3SStarPlus; + n_rk_stages = 4; + num_delta = 5; + } + else if (rk_method_string == "RK4_3_5_3SStar"){ + runge_kutta_method = RKMethodEnum::RK4_3_5_3SStar; + n_rk_stages = 5; + num_delta = 7; + } + prm.enter_subsection("rrk root solver"); { const std::string output_string_rrk = prm.get("rrk_root_solver_output"); diff --git a/src/parameters/parameters_ode_solver.h b/src/parameters/parameters_ode_solver.h index cf40521f7..bab49379e 100644 --- a/src/parameters/parameters_ode_solver.h +++ b/src/parameters/parameters_ode_solver.h @@ -62,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, + RK3_2_5F_3SStarPlus }; RKMethodEnum runge_kutta_method; ///< Runge-kutta method. int n_rk_stages; ///< Number of stages for an RK method; assigned based on runge_kutta_method int rk_order; ///< Order of the RK method; assigned based on runge_kutta_method + int num_delta; /// Flag to signal that automatic differentiation (AD) matrix dRdW must be allocated bool allocate_matrix_dRdW; From 40d4e00dcd9eb876f05dfc6e32ab350ac0e71003 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Wed, 3 Jul 2024 14:10:02 -0400 Subject: [PATCH 15/18] low storage algorithms, need to fix coefficients for RK3_2_5F_3SStarPlus --- src/flow_solver/flow_solver.cpp | 4 +- .../low_storage_runge_kutta_ode_solver.cpp | 58 +++++++++++++------ src/ode_solver/ode_solver_factory.cpp | 17 +++++- .../low_storage_rk_tableau_base.cpp | 3 +- .../low_storage_runge_kutta_methods.cpp | 18 ++++-- src/parameters/parameters_ode_solver.cpp | 10 +++- ...me_refinement_study_advection_explicit.prm | 4 +- 7 files changed, 81 insertions(+), 33 deletions(-) diff --git a/src/flow_solver/flow_solver.cpp b/src/flow_solver/flow_solver.cpp index 3ee9fb67b..426b62302 100644 --- a/src/flow_solver/flow_solver.cpp +++ b/src/flow_solver/flow_solver.cpp @@ -502,8 +502,8 @@ int FlowSolver::run() const // update time step in flow_solver_case 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; + //next_time_step = ode_solver->err_time_step(time_step,false); + // pcout << "next time step " << next_time_step; // advance solution // ode_solver->err_time_step(time_step); // pseudotime==false diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 38b5dd1de..6779c3067 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -26,51 +26,67 @@ void LowStorageRungeKuttaODESolver::step_in_time dealii::LinearAlgebra::distributed::Vector storage_register_2; storage_register_2.reinit(this->solution_update); storage_register_2*=0; - dealii::LinearAlgebra::distributed::Vector storage_register_3; - storage_register_3.reinit(this->solution_update); - storage_register_3 = this->dg->solution; - dealii::LinearAlgebra::distributed::Vector storage_register_1 = storage_register_3; + dealii::LinearAlgebra::distributed::Vector storage_register_1; + storage_register_1.reinit(this->solution_update); + storage_register_1 = this->dg->solution; + dealii::LinearAlgebra::distributed::Vector storage_register_3 = storage_register_1; dealii::LinearAlgebra::distributed::Vector rhs = storage_register_1; dealii::LinearAlgebra::distributed::Vector storage_register_4; - - if (this->butcher_tableau->get_b_hat(n_rk_stages+1) != 0){ - storage_register_4.reinit(this->solution_update); - storage_register_4 = this->dg->solution; + //this->pcout << this->butcher_tableau->get_b_hat(0) ; + //this->pcout << this->butcher_tableau->get_b_hat(1) ; + //this->pcout << this->butcher_tableau->get_b_hat(2) ; + //this->pcout << this->butcher_tableau->get_b_hat(3) ; + if (this->butcher_tableau->get_b_hat(0) != 0.0){ + this->pcout << "Initializing S4" << std::endl; + //storage_register_4.reinit(this->solution_update); + //storage_register_4 = this->dg->solution; + storage_register_4 = storage_register_1; } // int m = 5; 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; + //this->pcout << "Starting stage loop with i = " << i << std::endl; + //this->pcout << "deltai = " << this->butcher_tableau->get_delta(i-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); + //this->pcout << " gam1i = " << this->butcher_tableau->get_gamma(i, 0); storage_register_1.add(this->butcher_tableau->get_gamma(i, 1), storage_register_2); + //this->pcout << " gam2i = " << this->butcher_tableau->get_gamma(i, 1); storage_register_1.add(this->butcher_tableau->get_gamma(i, 2), storage_register_3); + //this->pcout << " gam3i = " << this->butcher_tableau->get_gamma(i, 2); rhs *= dt; storage_register_1.add(this->butcher_tableau->get_beta(i), rhs); - if (this->butcher_tableau->get_b_hat(n_rk_stages + 1) != 0){ - storage_register_4.add(this->butcher_tableau->get_b_hat(i), rhs); - } + this->pcout << " betai = " << this->butcher_tableau->get_beta(i); + 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 << "bhati = " << this->butcher_tableau->get_b_hat(i-1); + } + //this->pcout << std::endl; + // Check bhat (i) + // rhs = dt * f(S1) rhs = storage_register_1; - } - + // std::abort(); + //rhs = storage_register_1; // 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 < n_rk_stages +2; i++){ + for (int i = 0; i < n_rk_stages +2; 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); @@ -79,12 +95,13 @@ void LowStorageRungeKuttaODESolver::step_in_time // u_hat = s2 } else if (this->butcher_tableau->get_b_hat(n_rk_stages+1) != 0){ - storage_register_1 *= dt; - storage_register_4.add(this->butcher_tableau->get_b_hat(n_rk_stages + 1), storage_register_1); - // u_hat = s1 + rhs *= dt; + storage_register_4.add(this->butcher_tableau->get_b_hat(n_rk_stages + 1), rhs); + } 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){ @@ -110,7 +127,7 @@ void LowStorageRungeKuttaODESolver::step_in_time this->pcout << std::endl; // std::cout << "w " << w; - +*/ ++(this->current_iteration); this->current_time += dt; this->pcout << " Time is: " << this->current_time < >; template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; +template class LowStorageRungeKuttaODESolver >; #if PHILIP_DIM != 1 template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; template class LowStorageRungeKuttaODESolver >; + template class LowStorageRungeKuttaODESolver >; #endif } // ODESolver namespace diff --git a/src/ode_solver/ode_solver_factory.cpp b/src/ode_solver/ode_solver_factory.cpp index f2d6b77bf..50848a071 100644 --- a/src/ode_solver/ode_solver_factory.cpp +++ b/src/ode_solver/ode_solver_factory.cpp @@ -96,6 +96,7 @@ void ODESolverFactory::display_error_ode_solver_factory(Param else if (ode_solver_type == ODEEnum::rrk_explicit_solver) solver_string = "rrk_explicit"; else if (ode_solver_type == ODEEnum::pod_galerkin_solver) solver_string = "pod_galerkin"; else if (ode_solver_type == ODEEnum::pod_petrov_galerkin_solver) solver_string = "pod_petrov_galerkin"; + else if (ode_solver_type == ODEEnum::low_storage_runge_kutta_solver) solver_string = "low_storage_runge_kutta_solver"; else solver_string = "undefined"; dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0); @@ -122,14 +123,15 @@ template std::shared_ptr> ODESolverFactory::create_RungeKuttaODESolver(std::shared_ptr< DGBase > dg_input) { dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0); - std::shared_ptr> rk_tableau = create_RKTableau(dg_input); std::shared_ptr> RRK_object = create_RRKObject(dg_input, rk_tableau); + const int n_rk_stages = dg_input->all_parameters->ode_solver_param.n_rk_stages; using ODEEnum = Parameters::ODESolverParam::ODESolverEnum; const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type; if (ode_solver_type == ODEEnum::runge_kutta_solver || ode_solver_type == ODEEnum::rrk_explicit_solver) { + // Hard-coded templating of n_rk_stages because it is not known at compile time pcout << "Creating Runge Kutta ODE Solver with " << n_rk_stages << " stage(s)..." << std::endl; @@ -169,6 +171,9 @@ std::shared_ptr> ODESolverFactory>(dg_input,ls_rk_tableau,RRK_object); } + else if (n_rk_stages == 5){ + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + } else{ pcout << "Error: invalid number of stages. Aborting..." << std::endl; std::abort(); @@ -194,7 +199,7 @@ std::shared_ptr> ODESolverFactory> (n_rk_stages, num_delta, "RK4_3_5_3SStar"); //return std::make_shared> ("LSRK"); else { - pcout << "Error: invalid RK method. Aborting..." << std::endl; + pcout << "Error: invalid LS RK method. Aborting..." << std::endl; std::abort(); return nullptr; } @@ -226,6 +231,14 @@ std::shared_ptr> ODESolverFactory> (n_rk_stages, "2nd order diagonally-implicit (implicit)"); if (rk_method == RKMethodEnum::dirk_3_im) return std::make_shared> (n_rk_stages, "3nd order diagonally-implicit (implicit)"); else { + if (rk_method == RKMethodEnum::RK3_2_5F_3SStarPlus){ + pcout << "WARNING: returning dummy RK method!" << std::endl; + return std::make_shared> (n_rk_stages, "3rd order SSP (explicit)"); + } + if (rk_method == RKMethodEnum::RK4_3_5_3SStar) { + pcout << "WARNING: returning dummy RK method!" << std::endl; + return std::make_shared> (n_rk_stages, "3rd order SSP (explicit)"); + } pcout << "Error: invalid RK method. Aborting..." << std::endl; std::abort(); return nullptr; diff --git a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp index b358e04a0..4e37a20a1 100644 --- a/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp +++ b/src/ode_solver/runge_kutta_methods/low_storage_rk_tableau_base.cpp @@ -16,7 +16,7 @@ LowStorageRKTableauBase :: LowStorageRKTableauBase (const in this->butcher_tableau_gamma.reinit(n_rk_stages+1,3); this->butcher_tableau_beta.reinit(n_rk_stages+1); this->butcher_tableau_delta.reinit(num_delta); - this->butcher_tableau_b_hat.reinit(n_rk_stages+2); + this->butcher_tableau_b_hat.reinit(n_rk_stages+1); } template @@ -25,6 +25,7 @@ void LowStorageRKTableauBase :: set_tableau () set_gamma(); set_beta(); set_delta(); + set_b_hat(); pcout << "Assigned RK method: " << rk_method_string << std::endl; } diff --git a/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp b/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp index 56424b1b6..dcf58390f 100644 --- a/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp +++ b/src/ode_solver/runge_kutta_methods/low_storage_runge_kutta_methods.cpp @@ -10,7 +10,12 @@ namespace ODE { template void RK4_3_5_3SStar :: set_gamma() { - const double gamma[6][3] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {-0.497531095840104, 1.384996869124138, 0.0}, {1.010070514199942, 3.878155713328178, 0.0}, {-3.196559004608766,-2.324512951813145, 1.642598936063715}, {1.717835630267259, -0.514633322274467, 0.188295940828347}}; + const double gamma[6][3] = {{0.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {-0.497531095840104, 1.384996869124138, 0.0}, + {1.010070514199942, 3.878155713328178, 0.0}, + {-3.196559004608766,-2.324512951813145, 1.642598936063715}, + {1.717835630267259, -0.514633322274467, 0.188295940828347}}; for (int i = 0; i < 6; i++){ for (int j = 0; j < 3; j++){ this->butcher_tableau_gamma[i][j] = gamma[i][j]; @@ -47,8 +52,13 @@ void RK4_3_5_3SStar :: set_b_hat() template void RK3_2_5F_3SStarPlus :: set_gamma() { - const double gamma[5][3] = {{0.0, 1.0, 0.0}, {0.2587771979725733308135192812685323706, 0.5528354909301389892439698870483746541, 0.0}, {-0.1324380360140723382965420909764953437, 0.6731871608203061824849561782794643600, 0.0}, {0.05056033948190826045833606441415585735, 0.2803103963297672407841316576323901761, 0.2752563273304676380891217287572780582}, {0.5670532000739313812633197158607642990, 0.5521525447020610386070346724931300367, -0.8950526174674033822276061734289327568}}; - for (int i = 0; i < 5; i++){ + const double gamma[6][3] = {{0.0, 0.0, 0.0}, // Ignored + {0.0, 1.0, 0.0}, // first loop + {0.2587771979725733308135192812685323706, 0.5528354909301389892439698870483746541, 0.0}, + {-0.1324380360140723382965420909764953437, 0.6731871608203061824849561782794643600, 0.0}, + {0.05056033948190826045833606441415585735, 0.2803103963297672407841316576323901761, 0.2752563273304676380891217287572780582}, + {0.5670532000739313812633197158607642990, 0.5521525447020610386070346724931300367, -0.8950526174674033822276061734289327568}}; + for (int i = 0; i < 6; i++){ for (int j = 0; j < 3; j++){ this->butcher_tableau_gamma[i][j] = gamma[i][j]; } @@ -58,7 +68,7 @@ void RK3_2_5F_3SStarPlus :: set_gamma() template void RK3_2_5F_3SStarPlus :: set_beta() { - const double beta[5] = {0.11479359710235412, 0.089334428531133159, 0.43558710250086169, 0.24735761882014512, 0.058819144221557401}; + const double beta[6] = {0.0, 0.11479359710235412, 0.089334428531133159, 0.43558710250086169, 0.24735761882014512, 0.11292725304550591}; this->butcher_tableau_beta.fill(beta); } diff --git a/src/parameters/parameters_ode_solver.cpp b/src/parameters/parameters_ode_solver.cpp index 76b2fd23f..b2707efec 100644 --- a/src/parameters/parameters_ode_solver.cpp +++ b/src/parameters/parameters_ode_solver.cpp @@ -104,7 +104,9 @@ void ODESolverParam::declare_parameters (dealii::ParameterHandler &prm) " euler_ex | " " euler_im | " " dirk_2_im | " - " 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." "Choices are " " ."); + " dirk_3_im | " + " RK4_3_5_3SStar" + " RK3_2_5F_3SStarPlus | >."); prm.enter_subsection("rrk root solver"); { prm.declare_entry("rrk_root_solver_output", "quiet", @@ -215,7 +219,7 @@ void ODESolverParam::parse_parameters (dealii::ParameterHandler &prm) } else if (rk_method_string == "RK3_2_5F_3SStarPlus"){ runge_kutta_method = RKMethodEnum::RK3_2_5F_3SStarPlus; - n_rk_stages = 4; + n_rk_stages = 5; num_delta = 5; } else if (rk_method_string == "RK4_3_5_3SStar"){ diff --git a/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm index 9dea8271a..93c5c025b 100644 --- a/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm +++ b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm @@ -11,10 +11,10 @@ set use_periodic_bc = true # ODE solver subsection ODE solver - set ode_solver_type = runge_kutta + set ode_solver_type = low_storage_runge_kutta set output_solution_every_dt_time_intervals = 0.1 set initial_time_step = 2.5E-3 - set runge_kutta_method = ssprk3_ex + set runge_kutta_method = RK3_2_5F_3SStarPlus end subsection manufactured solution convergence study From 659719bb6e68026876769c7f09955c853894e908 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Wed, 3 Jul 2024 15:39:24 -0400 Subject: [PATCH 16/18] correct convergence for both algorithms (3S* and 3S*+) --- src/flow_solver/flow_solver.cpp | 2 +- .../low_storage_runge_kutta_ode_solver.cpp | 15 ++++++++------- .../low_storage_runge_kutta_methods.cpp | 9 +++++++++ .../euler_inviscid_TGV_largest_stable_CFL.prm | 1 + .../time_refinement_study_advection_explicit.prm | 5 +++-- 5 files changed, 22 insertions(+), 10 deletions(-) diff --git a/src/flow_solver/flow_solver.cpp b/src/flow_solver/flow_solver.cpp index 426b62302..f65b3afc0 100644 --- a/src/flow_solver/flow_solver.cpp +++ b/src/flow_solver/flow_solver.cpp @@ -502,7 +502,7 @@ int FlowSolver::run() const // update time step in flow_solver_case ode_solver->step_in_time(time_step,false); pcout << "time step " << time_step; - //next_time_step = ode_solver->err_time_step(time_step,false); + next_time_step = ode_solver->err_time_step(time_step,false); // pcout << "next time step " << next_time_step; // advance solution // ode_solver->err_time_step(time_step); // pseudotime==false diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 6779c3067..99cf3d493 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -49,12 +49,12 @@ void LowStorageRungeKuttaODESolver::step_in_time // int m = 5; double sum_delta = 0; - /* - double atol = 0.001; - double rtol = 0.001; + + double atol = 0.006; + double rtol = 0.006; 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; //this->pcout << "Starting stage loop with i = " << i << std::endl; @@ -83,7 +83,6 @@ void LowStorageRungeKuttaODESolver::step_in_time rhs = storage_register_1; } // std::abort(); - //rhs = storage_register_1; // 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 < n_rk_stages +2; i++){ //change n_rk_stages+2 to num_delta on this line @@ -92,6 +91,7 @@ void LowStorageRungeKuttaODESolver::step_in_time 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; + //this->pcout << "s2 " << sum_delta; // u_hat = s2 } else if (this->butcher_tableau->get_b_hat(n_rk_stages+1) != 0){ @@ -101,7 +101,7 @@ void LowStorageRungeKuttaODESolver::step_in_time } 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){ @@ -109,6 +109,7 @@ void LowStorageRungeKuttaODESolver::step_in_time for (dealii::LinearAlgebra::distributed::Vector::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); + //this->pcout << "w " << w; } } else if (this->butcher_tableau->get_b_hat(1) != 0){ @@ -127,7 +128,7 @@ void LowStorageRungeKuttaODESolver::step_in_time this->pcout << std::endl; // std::cout << "w " << w; -*/ + ++(this->current_iteration); this->current_time += dt; this->pcout << " Time is: " << this->current_time < :: set_gamma() } } +/* template void RK3_2_5F_3SStarPlus :: set_beta() { const double beta[6] = {0.0, 0.11479359710235412, 0.089334428531133159, 0.43558710250086169, 0.24735761882014512, 0.11292725304550591}; this->butcher_tableau_beta.fill(beta); } +*/ + +template +void RK3_2_5F_3SStarPlus :: set_beta() +{ + const double beta[6] = {0.0, 0.2300298624518076223899418286314123354, 0.3021434166948288809034402119555380003, 0.8025606185416310937583009085873554681, 0.4362158943603440930655148245148766471, 0.11292725304550591}; + this->butcher_tableau_beta.fill(beta); +} template void RK3_2_5F_3SStarPlus :: set_delta() diff --git a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm index dc64e5491..d79ea963d 100644 --- a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm +++ b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm @@ -27,6 +27,7 @@ subsection ODE solver set ode_solver_type = low_storage_runge_kutta #set ode_solver_type = rrk_explicit #set runge_kutta_method = ssprk3_ex + set runge_kutta_method = RK4_3_5_3SStar subsection rrk root solver set rrk_root_solver_output = verbose end diff --git a/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm index 93c5c025b..8f33cd44c 100644 --- a/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm +++ b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm @@ -13,7 +13,8 @@ set use_periodic_bc = true subsection ODE solver set ode_solver_type = low_storage_runge_kutta set output_solution_every_dt_time_intervals = 0.1 - set initial_time_step = 2.5E-3 + set initial_time_step = 5E-3 + #set runge_kutta_method = RK4_3_5_3SStar set runge_kutta_method = RK3_2_5F_3SStarPlus end @@ -26,7 +27,7 @@ end subsection time_refinement_study set number_of_times_to_solve = 4 - set refinement_ratio = 0.5 + set refinement_ratio = 0.8 end subsection flow_solver From 4a99eac51ccc97ec2bddb3533d8ab6f17c9b9abe Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Fri, 5 Jul 2024 16:14:26 -0400 Subject: [PATCH 17/18] add tolerances as parameters --- .../flow_solver_cases/periodic_turbulence.cpp | 15 +++++- .../low_storage_runge_kutta_ode_solver.cpp | 50 +++++++------------ .../low_storage_runge_kutta_ode_solver.h | 5 ++ src/ode_solver/ode_solver_factory.cpp | 14 +++--- src/parameters/parameters_flow_solver.cpp | 1 - src/parameters/parameters_flow_solver.h | 2 + src/parameters/parameters_linear_solver.cpp | 4 ++ src/parameters/parameters_linear_solver.h | 4 +- src/parameters/parameters_ode_solver.cpp | 2 + ...nviscid_isentropic_vortex_h_refinement.prm | 5 +- .../euler_inviscid_TGV_largest_stable_CFL.prm | 2 +- ...me_refinement_study_advection_explicit.prm | 6 ++- 12 files changed, 64 insertions(+), 46 deletions(-) diff --git a/src/flow_solver/flow_solver_cases/periodic_turbulence.cpp b/src/flow_solver/flow_solver_cases/periodic_turbulence.cpp index 50fd7b120..7e00d8d32 100644 --- a/src/flow_solver/flow_solver_cases/periodic_turbulence.cpp +++ b/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 #include #include @@ -686,7 +686,15 @@ void PeriodicTurbulence::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_solver_low_storage = std::dynamic_pointer_cast>(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::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); diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp index 99cf3d493..12710a3d4 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.cpp @@ -17,7 +17,7 @@ LowStorageRungeKuttaODESolver::LowStorageRungeKu } template -void LowStorageRungeKuttaODESolver::step_in_time (real dt, const bool pseudotime) +void LowStorageRungeKuttaODESolver::step_in_time (real dt, const bool pseudotime) { this->original_time_step = dt; this->solution_update = this->dg->solution; //storing u_n @@ -34,48 +34,32 @@ void LowStorageRungeKuttaODESolver::step_in_time dealii::LinearAlgebra::distributed::Vector rhs = storage_register_1; dealii::LinearAlgebra::distributed::Vector storage_register_4; - //this->pcout << this->butcher_tableau->get_b_hat(0) ; - //this->pcout << this->butcher_tableau->get_b_hat(1) ; - //this->pcout << this->butcher_tableau->get_b_hat(2) ; - //this->pcout << this->butcher_tableau->get_b_hat(3) ; + if (this->butcher_tableau->get_b_hat(0) != 0.0){ - this->pcout << "Initializing S4" << std::endl; - //storage_register_4.reinit(this->solution_update); - //storage_register_4 = this->dg->solution; storage_register_4 = storage_register_1; } - - - // int m = 5; double sum_delta = 0; - double atol = 0.006; - double rtol = 0.006; + //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; - //this->pcout << "Starting stage loop with i = " << i << std::endl; - //this->pcout << "deltai = " << this->butcher_tableau->get_delta(i-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); - //this->pcout << " gam1i = " << this->butcher_tableau->get_gamma(i, 0); storage_register_1.add(this->butcher_tableau->get_gamma(i, 1), storage_register_2); - //this->pcout << " gam2i = " << this->butcher_tableau->get_gamma(i, 1); storage_register_1.add(this->butcher_tableau->get_gamma(i, 2), storage_register_3); - //this->pcout << " gam3i = " << this->butcher_tableau->get_gamma(i, 2); rhs *= dt; storage_register_1.add(this->butcher_tableau->get_beta(i), rhs); - this->pcout << " betai = " << this->butcher_tableau->get_beta(i); 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 << "bhati = " << this->butcher_tableau->get_b_hat(i-1); } //this->pcout << std::endl; // Check bhat (i) @@ -85,23 +69,28 @@ void LowStorageRungeKuttaODESolver::step_in_time // 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 < n_rk_stages +2; i++){ //change n_rk_stages+2 to num_delta on this line + 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; - //this->pcout << "s2 " << sum_delta; // u_hat = s2 } - else if (this->butcher_tableau->get_b_hat(n_rk_stages+1) != 0){ + + + // 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 + 1), rhs); + 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){ @@ -109,7 +98,6 @@ void LowStorageRungeKuttaODESolver::step_in_time for (dealii::LinearAlgebra::distributed::Vector::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); - //this->pcout << "w " << w; } } else if (this->butcher_tableau->get_b_hat(1) != 0){ @@ -140,17 +128,17 @@ template double LowStorageRungeKuttaODESolver::err_time_step (real dt, const bool pseudotime) { (void) pseudotime; - int q_hat = 3; - int k = q_hat +1; + //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]/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; @@ -158,7 +146,7 @@ double LowStorageRungeKuttaODESolver::err_time_s template -void LowStorageRungeKuttaODESolver::allocate_ode_system () +void LowStorageRungeKuttaODESolver::allocate_ode_system () { this->pcout << "Allocating ODE system..." << std::flush; this->solution_update.reinit(this->dg->right_hand_side); diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.h b/src/ode_solver/low_storage_runge_kutta_ode_solver.h index 83b210e6c..eeae7a4c7 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.h +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.h @@ -54,6 +54,11 @@ class LowStorageRungeKuttaODESolver: public ODESolverBase real w; double epsilon[3]; + int num_delta; + int rk_order; + double atol; + double rtol; + //double beta_controller[3]; }; } // ODE namespace diff --git a/src/ode_solver/ode_solver_factory.cpp b/src/ode_solver/ode_solver_factory.cpp index 50848a071..bea4de2ce 100644 --- a/src/ode_solver/ode_solver_factory.cpp +++ b/src/ode_solver/ode_solver_factory.cpp @@ -128,6 +128,7 @@ std::shared_ptr> ODESolverFactoryall_parameters->ode_solver_param.n_rk_stages; + //const int num_delta = dg_input->all_parameters->ode_solver_param.num_delta; using ODEEnum = Parameters::ODESolverParam::ODESolverEnum; const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type; if (ode_solver_type == ODEEnum::runge_kutta_solver || ode_solver_type == ODEEnum::rrk_explicit_solver) { @@ -159,20 +160,19 @@ std::shared_ptr> ODESolverFactory>(dg_input,ls_rk_tableau,RRK_object); + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); } else if (n_rk_stages == 2){ - - return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); } else if (n_rk_stages == 3){ - return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); } - else if (n_rk_stages == 4){ - return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + else if (n_rk_stages == 4){ + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); } else if (n_rk_stages == 5){ - return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); + return std::make_shared>(dg_input,ls_rk_tableau,RRK_object); } else{ pcout << "Error: invalid number of stages. Aborting..." << std::endl; diff --git a/src/parameters/parameters_flow_solver.cpp b/src/parameters/parameters_flow_solver.cpp index bd2772633..97f8a878b 100644 --- a/src/parameters/parameters_flow_solver.cpp +++ b/src/parameters/parameters_flow_solver.cpp @@ -349,7 +349,6 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm) max_poly_degree_for_adaptation = prm.get_integer("max_poly_degree_for_adaptation"); // -- set value to poly_degree if it is the default value if(max_poly_degree_for_adaptation == 0) max_poly_degree_for_adaptation = poly_degree; - final_time = prm.get_double("final_time"); constant_time_step = prm.get_double("constant_time_step"); courant_friedrichs_lewy_number = prm.get_double("courant_friedrichs_lewy_number"); diff --git a/src/parameters/parameters_flow_solver.h b/src/parameters/parameters_flow_solver.h index 5fa89df7e..4f9a4e63e 100644 --- a/src/parameters/parameters_flow_solver.h +++ b/src/parameters/parameters_flow_solver.h @@ -39,6 +39,8 @@ class FlowSolverParam unsigned int poly_degree; ///< Polynomial order (P) of the basis functions for DG. unsigned int max_poly_degree_for_adaptation; ///< Maximum polynomial order of the DG basis functions for adaptation. + + double final_time; ///< Final solution time double constant_time_step; ///< Constant time step double courant_friedrichs_lewy_number; ///< Courant-Friedrich-Lewy (CFL) number for constant time step diff --git a/src/parameters/parameters_linear_solver.cpp b/src/parameters/parameters_linear_solver.cpp index 27c2da69d..91cd72813 100644 --- a/src/parameters/parameters_linear_solver.cpp +++ b/src/parameters/parameters_linear_solver.cpp @@ -82,6 +82,10 @@ void LinearSolverParam ::parse_parameters (dealii::ParameterHandler &prm) { prm.enter_subsection("linear solver"); { + atol = prm.get_double("atol"); + rtol = prm.get_double("rtol"); + //beta_controller = prm.get_double("beta_controller"); + const std::string output_string = prm.get("linear_solver_output"); if (output_string == "verbose") linear_solver_output = verbose; if (output_string == "quiet") linear_solver_output = quiet; diff --git a/src/parameters/parameters_linear_solver.h b/src/parameters/parameters_linear_solver.h index 43629be6e..a24d3f774 100644 --- a/src/parameters/parameters_linear_solver.h +++ b/src/parameters/parameters_linear_solver.h @@ -24,7 +24,9 @@ class LinearSolverParam */ OutputEnum linear_solver_output; ///< quiet or verbose. LinearSolverEnum linear_solver_type; ///< direct or gmres. - + double atol; + double rtol; + //double beta_controller[3]; // GMRES options double ilut_drop; ///< Threshold to drop terms close to zero. double ilut_rtol; ///< Multiplies diagonal by ilut_rtol for more diagonal dominance. diff --git a/src/parameters/parameters_ode_solver.cpp b/src/parameters/parameters_ode_solver.cpp index b2707efec..ad2ac7973 100644 --- a/src/parameters/parameters_ode_solver.cpp +++ b/src/parameters/parameters_ode_solver.cpp @@ -221,11 +221,13 @@ void ODESolverParam::parse_parameters (dealii::ParameterHandler &prm) runge_kutta_method = RKMethodEnum::RK3_2_5F_3SStarPlus; n_rk_stages = 5; num_delta = 5; + rk_order = 3; } else if (rk_method_string == "RK4_3_5_3SStar"){ runge_kutta_method = RKMethodEnum::RK4_3_5_3SStar; n_rk_stages = 5; num_delta = 7; + rk_order = 4; } prm.enter_subsection("rrk root solver"); diff --git a/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm b/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm index 129130330..ff684e774 100644 --- a/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm +++ b/tests/integration_tests_control_files/isentropic_vortex_integration/2D_inviscid_isentropic_vortex_h_refinement.prm @@ -26,10 +26,13 @@ set two_point_num_flux_type = IR subsection ODE solver set ode_output = quiet set ode_solver_type = low_storage_runge_kutta -# set runge_kutta_method = rk4_ex + set runge_kutta_method = RK3_2_5F_3SStarPlus end subsection linear solver + set atol = 0.001 + set rtol = 0.001 + #set beta_controller = {0.70, -0.23, 0} set linear_solver_output = quiet end diff --git a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm index d79ea963d..9c4bef708 100644 --- a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm +++ b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm @@ -27,7 +27,7 @@ subsection ODE solver set ode_solver_type = low_storage_runge_kutta #set ode_solver_type = rrk_explicit #set runge_kutta_method = ssprk3_ex - set runge_kutta_method = RK4_3_5_3SStar + set runge_kutta_method = RK3_2_5F_3SStarPlus subsection rrk root solver set rrk_root_solver_output = verbose end diff --git a/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm index 8f33cd44c..24e1180b4 100644 --- a/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm +++ b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm @@ -13,7 +13,8 @@ set use_periodic_bc = true subsection ODE solver set ode_solver_type = low_storage_runge_kutta set output_solution_every_dt_time_intervals = 0.1 - set initial_time_step = 5E-3 + #set initial_time_step = 5E-3 (convergence test) + set initial_time_step = 2.5E-3 #set runge_kutta_method = RK4_3_5_3SStar set runge_kutta_method = RK3_2_5F_3SStarPlus end @@ -26,7 +27,8 @@ end subsection time_refinement_study - set number_of_times_to_solve = 4 + #set number_of_times_to_solve = 4 + set number_of_times_to_solve = 1 set refinement_ratio = 0.8 end From ac2b44672fb4431509c2bef0ba5627f117ef5d30 Mon Sep 17 00:00:00 2001 From: jocelynpoon Date: Mon, 8 Jul 2024 14:56:48 -0400 Subject: [PATCH 18/18] tolerance parameters work --- src/ode_solver/low_storage_runge_kutta_ode_solver.h | 9 ++++++++- src/parameters/parameters_linear_solver.cpp | 13 +++++++++++++ src/parameters/parameters_linear_solver.h | 6 +++--- .../euler_inviscid_TGV_largest_stable_CFL.prm | 5 +++++ 4 files changed, 29 insertions(+), 4 deletions(-) diff --git a/src/ode_solver/low_storage_runge_kutta_ode_solver.h b/src/ode_solver/low_storage_runge_kutta_ode_solver.h index eeae7a4c7..d9c12df14 100644 --- a/src/ode_solver/low_storage_runge_kutta_ode_solver.h +++ b/src/ode_solver/low_storage_runge_kutta_ode_solver.h @@ -50,7 +50,14 @@ class LowStorageRungeKuttaODESolver: public ODESolverBase /// Indicator for zero diagonal elements; used to toggle implicit solve. std::vector butcher_tableau_aii_is_zero; - +/* + std::vector> storage_register_2; + std::vector> storage_register_1; + std::vector> storage_register_3; + std::vector> storage_register_4; + std::vector> rhs; +*/ + real w; double epsilon[3]; diff --git a/src/parameters/parameters_linear_solver.cpp b/src/parameters/parameters_linear_solver.cpp index 91cd72813..fbfdf6f1e 100644 --- a/src/parameters/parameters_linear_solver.cpp +++ b/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", + dealii::Patterns::Double(), + "Absolute Tolerance"); + + prm.declare_entry("rtol", "0.001", + dealii::Patterns::Double(), + "Relative Tolerance"); +/* + prm.declare_entry("beta_controller", "{0.70, -0.23, 0}", + dealii::Patterns::Double(), + "Beta controller"); +*/ prm.declare_entry("linear_solver_output", "quiet", dealii::Patterns::Selection("quiet|verbose"), "State whether output from linear solver should be printed. " @@ -84,6 +96,7 @@ void LinearSolverParam ::parse_parameters (dealii::ParameterHandler &prm) { atol = prm.get_double("atol"); rtol = prm.get_double("rtol"); + //beta_controller = prm.get_double("beta_controller"); const std::string output_string = prm.get("linear_solver_output"); diff --git a/src/parameters/parameters_linear_solver.h b/src/parameters/parameters_linear_solver.h index a24d3f774..872a1a67b 100644 --- a/src/parameters/parameters_linear_solver.h +++ b/src/parameters/parameters_linear_solver.h @@ -24,9 +24,9 @@ class LinearSolverParam */ OutputEnum linear_solver_output; ///< quiet or verbose. LinearSolverEnum linear_solver_type; ///< direct or gmres. - double atol; - double rtol; - //double beta_controller[3]; + double atol; // Absolute tolerance + double rtol; // Relative tolerance + double beta_controller[3]; // GMRES options double ilut_drop; ///< Threshold to drop terms close to zero. double ilut_rtol; ///< Multiplies diagonal by ilut_rtol for more diagonal dominance. diff --git a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm index 9c4bef708..bc6d27333 100644 --- a/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm +++ b/tests/integration_tests_control_files/relaxation_runge_kutta_checks/euler_inviscid_TGV_largest_stable_CFL.prm @@ -21,6 +21,11 @@ set use_periodic_bc = true set conv_num_flux = two_point_flux set two_point_num_flux_type = Ra +subsection linear solver + set atol = 0.001 + set rtol = 0.001 +end + # ODE solver subsection ODE solver set ode_output = quiet