1 #include "ode_solver_factory.h"     2 #include "parameters/all_parameters.h"     3 #include "ode_solver_base.h"     4 #include "runge_kutta_ode_solver.h"     5 #include "low_storage_runge_kutta_ode_solver.h"     6 #include "implicit_ode_solver.h"     7 #include "relaxation_runge_kutta/algebraic_rrk_ode_solver.h"     8 #include "relaxation_runge_kutta/root_finding_rrk_ode_solver.h"     9 #include "pod_galerkin_ode_solver.h"    10 #include "pod_petrov_galerkin_ode_solver.h"    11 #include "pod_galerkin_runge_kutta_ode_solver.h"    12 #include "hyper_reduced_petrov_galerkin_ode_solver.h"    13 #include <deal.II/distributed/solution_transfer.h>    14 #include "runge_kutta_methods/runge_kutta_methods.h"    15 #include "runge_kutta_methods/low_storage_runge_kutta_methods.h"    16 #include "relaxation_runge_kutta/empty_RRK_base.h"    21 template <
int dim, 
typename real, 
typename MeshType>
    24     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
    25     pcout << 
"Creating ODE Solver..." << std::endl;
    27     const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
    28     if((ode_solver_type == ODEEnum::runge_kutta_solver)||(ode_solver_type == ODEEnum::rrk_explicit_solver || ode_solver_type == ODEEnum::low_storage_runge_kutta_solver))     
    29         return create_RungeKuttaODESolver(dg_input); 
    30     if(ode_solver_type == ODEEnum::implicit_solver)         
    31         return std::make_shared<ImplicitODESolver<dim,real,MeshType>>(dg_input);
    33         display_error_ode_solver_factory(ode_solver_type, 
false);
    38 template <
int dim, 
typename real, 
typename MeshType>
    41     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
    42     pcout << 
"Creating ODE Solver..." << std::endl;
    44     const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
    45     if(ode_solver_type == ODEEnum::pod_galerkin_solver)
    46         return std::make_shared<PODGalerkinODESolver<dim,real,MeshType>>(dg_input, pod);
    47     if(ode_solver_type == ODEEnum::pod_petrov_galerkin_solver) 
    48         return std::make_shared<PODPetrovGalerkinODESolver<dim,real,MeshType>>(dg_input, pod);
    49     if(ode_solver_type == ODEEnum::pod_galerkin_runge_kutta_solver)
    50         return create_RungeKuttaODESolver(dg_input, pod);
    52         display_error_ode_solver_factory(ode_solver_type, 
true);
    57 template <
int dim, 
typename real, 
typename MeshType>
    60     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
    61     pcout << 
"Creating ODE Solver..." << std::endl;
    63     const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
    64     if(ode_solver_type == ODEEnum::hyper_reduced_petrov_galerkin_solver) 
    65         return std::make_shared<HyperReducedODESolver<dim,real,MeshType>>(dg_input, pod, weights);
    67         display_error_ode_solver_factory(ode_solver_type, 
true);
    72 template <
int dim, 
typename real, 
typename MeshType>
    75     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
    76     pcout << 
"Creating ODE Solver..." << std::endl;
    78     if((ode_solver_type == ODEEnum::runge_kutta_solver)||(ode_solver_type == ODEEnum::rrk_explicit_solver))     
    79         return create_RungeKuttaODESolver(dg_input);
    80     if(ode_solver_type == ODEEnum::implicit_solver)         
    81         return std::make_shared<ImplicitODESolver<dim,real,MeshType>>(dg_input);
    83         display_error_ode_solver_factory(ode_solver_type, 
false);
    88 template <
int dim, 
typename real, 
typename MeshType>
    91     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
    92     pcout << 
"Creating ODE Solver..." << std::endl;
    94     if(ode_solver_type == ODEEnum::pod_galerkin_solver) 
    95         return std::make_shared<PODGalerkinODESolver<dim,real,MeshType>>(dg_input, pod);
    96     if(ode_solver_type == ODEEnum::pod_petrov_galerkin_solver) 
    97         return std::make_shared<PODPetrovGalerkinODESolver<dim,real,MeshType>>(dg_input, pod);
    98     if(ode_solver_type == ODEEnum::pod_galerkin_runge_kutta_solver)
    99         return create_RungeKuttaODESolver(dg_input, pod);
   101         display_error_ode_solver_factory(ode_solver_type, 
true);
   106 template <
int dim, 
typename real, 
typename MeshType>
   109     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
   110     pcout << 
"Creating ODE Solver..." << std::endl;
   112     if(ode_solver_type == ODEEnum::hyper_reduced_petrov_galerkin_solver) 
   113         return std::make_shared<HyperReducedODESolver<dim,real,MeshType>>(dg_input, pod, weights);
   115         display_error_ode_solver_factory(ode_solver_type, 
true);
   121 template <
int dim, 
typename real, 
typename MeshType>
   126     std::string solver_string;    
   127     if (ode_solver_type == ODEEnum::runge_kutta_solver)                 solver_string = 
"runge_kutta";
   128     else if (ode_solver_type == ODEEnum::implicit_solver)               solver_string = 
"implicit";
   129     else if (ode_solver_type == ODEEnum::rrk_explicit_solver)           solver_string = 
"rrk_explicit";
   130     else if (ode_solver_type == ODEEnum::pod_galerkin_solver)           solver_string = 
"pod_galerkin";
   131     else if (ode_solver_type == ODEEnum::pod_petrov_galerkin_solver)    solver_string = 
"pod_petrov_galerkin";
   132     else if (ode_solver_type == ODEEnum::hyper_reduced_petrov_galerkin_solver)    
   133                                                                         solver_string = 
"hyper_reduced_petrov_galerkin";
   134     else if (ode_solver_type == ODEEnum::low_storage_runge_kutta_solver)    
   135                                                                         solver_string = 
"low_storage_runge_kutta_solver";
   136     else if (ode_solver_type == ODEEnum::pod_galerkin_runge_kutta_solver)
   137                                                                         solver_string = 
"pod_galerkin_runge_kutta";
   138     else solver_string = 
"undefined";
   141     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
   142     pcout << 
"********************************************************************" << std::endl;
   143     pcout << 
"Can't create ODE solver since solver type is not clear." << std::endl;
   144     pcout << 
"Solver type specified: " << solver_string << std::endl;
   145     pcout << 
"Solver type possible: " << std::endl;
   147         pcout <<  
"pod_galerkin" << std::endl;
   148         pcout <<  
"pod_petrov_galerkin" << std::endl;
   149         pcout <<  
"pod_galerkin_runge_kutta" << std::endl;
   152         pcout <<  
"runge_kutta" << std::endl;
   153         pcout <<  
"implicit" << std::endl;
   154         pcout <<  
"rrk_explicit" << std::endl;
   155         pcout << 
"    With rrk_explicit only being valid for " <<std::endl;
   156         pcout << 
"    pde_type = burgers, flux_nodes_type = GLL, overintegration = 0, and dim = 1" <<std::endl;
   158     pcout << 
"********************************************************************" << std::endl;
   162 template <
int dim, 
typename real, 
typename MeshType>
   165     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
   166     std::shared_ptr<RKTableauBase<dim,real,MeshType>> rk_tableau = create_RKTableau(dg_input);
   167     std::shared_ptr<EmptyRRKBase<dim,real,MeshType>> RRK_object = create_RRKObject(dg_input, rk_tableau);
   170     const int n_rk_stages = dg_input->all_parameters->ode_solver_param.n_rk_stages;
   172     const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
   173     if (ode_solver_type == ODEEnum::runge_kutta_solver || ode_solver_type == ODEEnum::rrk_explicit_solver) {
   178         pcout << 
"Creating Runge Kutta ODE Solver with "    179               << n_rk_stages << 
" stage(s)..." << std::endl;
   180         if (n_rk_stages == 1){
   181             return std::make_shared<RungeKuttaODESolver<dim,real,1,MeshType>>(dg_input,rk_tableau_butcher,RRK_object);
   183         else if (n_rk_stages == 2){
   184             return std::make_shared<RungeKuttaODESolver<dim,real,2,MeshType>>(dg_input,rk_tableau_butcher,RRK_object);
   186         else if (n_rk_stages == 3){
   187             return std::make_shared<RungeKuttaODESolver<dim,real,3,MeshType>>(dg_input,rk_tableau_butcher,RRK_object);
   189         else if (n_rk_stages == 4){
   190             return std::make_shared<RungeKuttaODESolver<dim,real,4,MeshType>>(dg_input,rk_tableau_butcher,RRK_object);
   193             pcout << 
"Error: invalid number of stages. Aborting..." << std::endl;
   197     } 
else if (ode_solver_type == ODEEnum::low_storage_runge_kutta_solver) {
   201         pcout << 
"Creating Low-Storage Runge Kutta ODE Solver with "    202               << n_rk_stages << 
" stage(s)..." << std::endl;
   203         if (n_rk_stages == 1){
   204             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,1, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   206         else if (n_rk_stages == 2){
   207             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,2, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   209         else if (n_rk_stages == 3){
   210             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,3, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   212         else if (n_rk_stages == 4){      
   213             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,4, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   215         else if (n_rk_stages == 5){
   216             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,5, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   218         else if (n_rk_stages == 9){
   219             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,9, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   221         else if (n_rk_stages == 10){
   222             return std::make_shared<LowStorageRungeKuttaODESolver<dim,real,10, MeshType>>(dg_input,ls_rk_tableau,RRK_object);
   225             pcout << 
"Error: invalid number of stages. Aborting..." << std::endl;
   230         display_error_ode_solver_factory(ode_solver_type, 
false);
   235 template <
int dim, 
typename real, 
typename MeshType> 
   238     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
   240     std::shared_ptr<RKTableauBase<dim,real,MeshType>> rk_tableau = create_RKTableau(dg_input);
   241     std::shared_ptr<EmptyRRKBase<dim,real,MeshType>> RRK_object = create_RRKObject(dg_input, rk_tableau);
   243     const int n_rk_stages = dg_input->all_parameters->ode_solver_param.n_rk_stages;
   245     const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
   246     if (ode_solver_type == ODEEnum::pod_galerkin_runge_kutta_solver) {
   250         if(dg_input->all_parameters->use_inverse_mass_on_the_fly){
   251             pcout   << 
"Not Implemented: use_inverse_mass_on_the_fly=true && ode_solver_type=pod_galerkin_rk_solver"   253                     << 
"Please set use_inverse_mass_on_the_fly=false and try again"   257         pcout << 
"Creating Galerkin Runge Kutta ODE Solver with "    258               << n_rk_stages << 
" stage(s)..." << std::endl;
   259         if (n_rk_stages == 1){
   260             return std::make_shared<PODGalerkinRungeKuttaODESolver<dim,real,1,MeshType>>(dg_input,rk_tableau_butcher,RRK_object,pod);
   262         else if (n_rk_stages == 2){
   263             return std::make_shared<PODGalerkinRungeKuttaODESolver<dim,real,2,MeshType>>(dg_input,rk_tableau_butcher,RRK_object,pod);
   265         else if (n_rk_stages == 3){
   266             return std::make_shared<PODGalerkinRungeKuttaODESolver<dim,real,3,MeshType>>(dg_input,rk_tableau_butcher,RRK_object,pod);
   268         else if (n_rk_stages == 4){
   269             return std::make_shared<PODGalerkinRungeKuttaODESolver<dim,real,4,MeshType>>(dg_input,rk_tableau_butcher,RRK_object,pod);
   272             pcout << 
"Error: invalid number of stages. Aborting..." << std::endl;
   278         display_error_ode_solver_factory(ode_solver_type, 
false);
   283 template <
int dim, 
typename real, 
typename MeshType>
   286     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
   288     const RKMethodEnum rk_method = dg_input->all_parameters->ode_solver_param.runge_kutta_method;
   290     const int n_rk_stages = dg_input->all_parameters->ode_solver_param.n_rk_stages;
   292     if (rk_method == RKMethodEnum::ssprk3_ex)   
return std::make_shared<SSPRK3Explicit<dim, real, MeshType>> (n_rk_stages, 
"3rd order SSP (explicit)");
   293     if (rk_method == RKMethodEnum::rk4_ex)      
return std::make_shared<RK4Explicit<dim, real, MeshType>>    (n_rk_stages, 
"4th order classical RK (explicit)");
   294     if (rk_method == RKMethodEnum::heun2_ex)      
return std::make_shared<HeunExplicit<dim, real, MeshType>>    (n_rk_stages, 
"2nd order Heun's method (explicit)");
   295     if (rk_method == RKMethodEnum::euler_ex) {
   297         ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
   298         if (ode_solver_type == ODEEnum::rrk_explicit_solver) {
   300             pcout << 
"Error: RRK is not valid for Forward Euler. Aborting..." << std::endl;
   303         } 
else return std::make_shared<EulerExplicit<dim, real, MeshType>>  (n_rk_stages, 
"Forward Euler (explicit)");
   305     if (rk_method == RKMethodEnum::euler_im)    
return std::make_shared<EulerImplicit<dim, real, MeshType>>  (n_rk_stages, 
"Implicit Euler (implicit)");
   306     if (rk_method == RKMethodEnum::dirk_2_im)   
return std::make_shared<DIRK2Implicit<dim, real, MeshType>>  (n_rk_stages, 
"2nd order diagonally-implicit (implicit)");
   307     if (rk_method == RKMethodEnum::dirk_3_im)   
return std::make_shared<DIRK3Implicit<dim, real, MeshType>>  (n_rk_stages, 
"3nd order diagonally-implicit (implicit)");
   310     const int num_delta = dg_input->all_parameters->ode_solver_param.num_delta;
   312     if (rk_method == RKMethodEnum::RK3_2_5F_3SStarPlus)         
return std::make_shared<RK3_2_5F_3SStarPlus<dim, real, MeshType>> (n_rk_stages, num_delta, 
"RK3_2_5F_3SStarPlus");
   313     if (rk_method == RKMethodEnum::RK4_3_5_3SStar)              
return std::make_shared<RK4_3_5_3SStar<dim, real, MeshType>>    (n_rk_stages, num_delta, 
"RK4_3_5_3SStar");
   314     if (rk_method == RKMethodEnum::RK4_3_9F_3SStarPlus)         
return std::make_shared<RK4_3_9F_3SStarPlus<dim, real, MeshType>>    (n_rk_stages, num_delta, 
"RK4_3_9F_3SStarPlus");
   315     if (rk_method == RKMethodEnum::RK5_4_10F_3SStarPlus)        
return std::make_shared<RK5_4_10F_3SStarPlus<dim, real, MeshType>>    (n_rk_stages, num_delta, 
"RK5_4_10F_3SStarPlus");
   317         pcout << 
"Error: invalid RK method. Aborting..." << std::endl;
   323 template <
int dim, 
typename real, 
typename MeshType>
   327     dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
   329     const ODEEnum ode_solver_type = dg_input->all_parameters->ode_solver_param.ode_solver_type;
   333     if ( (ode_solver_type == ODEEnum::runge_kutta_solver && dg_input->all_parameters->flow_solver_param.do_calculate_numerical_entropy)
   334             || ( !dg_input->all_parameters->ode_solver_param.use_relaxation_runge_kutta && dg_input->all_parameters->flow_solver_param.do_calculate_numerical_entropy )  ) {
   335             return std::make_shared<RKNumEntropy<dim,real,MeshType>>(rk_tableau_butcher);
   337     else if (dg_input->all_parameters->ode_solver_param.use_relaxation_runge_kutta){
   340         const PDEEnum pde_type = dg_input->all_parameters->pde_type;
   342         const NumFluxEnum two_point_num_flux_type = dg_input->all_parameters->two_point_num_flux_type;
   344         enum NumEntropyEnum {energy, nonlinear};
   345         NumEntropyEnum numerical_entropy_type;
   346         std::string rrk_type_string;
   347         if (pde_type == PDEEnum::burgers_inviscid){
   348             numerical_entropy_type = NumEntropyEnum::energy;
   349             rrk_type_string = 
"Algebraic";
   350         } 
else if ((pde_type == PDEEnum::euler || pde_type == PDEEnum::navier_stokes)
   351                     && (two_point_num_flux_type != NumFluxEnum::KG)){
   352             numerical_entropy_type = NumEntropyEnum::nonlinear;
   353             rrk_type_string = 
"Root-finding";
   355             pcout << 
"PDE type has no assigned numerical entropy variable. Aborting..." << std::endl;
   359         pcout << 
"Adding " << rrk_type_string << 
" Relaxation Runge Kutta to the ODE solver..." << std::endl;
   360         if (numerical_entropy_type==NumEntropyEnum::energy)
   361             return std::make_shared<AlgebraicRRKODESolver<dim,real,MeshType>>(rk_tableau_butcher);
   362         else if (numerical_entropy_type==NumEntropyEnum::nonlinear)
   363             return std::make_shared<RootFindingRRKODESolver<dim,real,MeshType>>(rk_tableau_butcher);
   366         return std::make_shared<EmptyRRKBase<dim,real,MeshType>> (rk_tableau_butcher);
 static std::shared_ptr< ODESolverBase< dim, real, MeshType > > create_ODESolver_manual(Parameters::ODESolverParam::ODESolverEnum ode_solver_type, std::shared_ptr< DGBase< dim, real, MeshType > > dg_input)
Creates either implicit or explicit ODE solver based on manual input (no POD basis given) ...
PartialDifferentialEquation
Possible Partial Differential Equations to solve. 
Files for the baseline physics. 
RKMethodEnum
Types of RK method (i.e., unique Butcher tableau) 
Base class for storing the RK method. 
ODESolverEnum
Types of ODE solver. 
Create specified ODE solver as ODESolverBase object. 
static std::shared_ptr< ODESolverBase< dim, real, MeshType > > create_RungeKuttaODESolver(std::shared_ptr< DGBase< dim, real, MeshType > > dg_input)
Creates an ODESolver object based on the specified RK method, including derived classes. 
TwoPointNumericalFlux
Two point numerical flux type for split form. 
static void display_error_ode_solver_factory(Parameters::ODESolverParam::ODESolverEnum ode_solver_type, bool reduced_order)
Output error message for Implicit and Explicit solver. 
static std::shared_ptr< RKTableauBase< dim, real, MeshType > > create_RKTableau(std::shared_ptr< DGBase< dim, real, MeshType > > dg_input)
Creates an RKTableau object based on the specified RK method. 
Base class for storing the RK method. 
DGBase is independent of the number of state variables. 
Base class for storing the RK method. 
static std::shared_ptr< ODESolverBase< dim, real, MeshType > > create_ODESolver(std::shared_ptr< DGBase< dim, real, MeshType > > dg_input)
Creates either implicit or explicit ODE solver based on parameter value(no POD basis given) ...
static std::shared_ptr< EmptyRRKBase< dim, real, MeshType > > create_RRKObject(std::shared_ptr< DGBase< dim, real, MeshType > > dg_input, std::shared_ptr< RKTableauBase< dim, real, MeshType >> rk_tableau)
Creates an RRK object with specified RRK type; if no RRK is being used, creates an RRK object with em...