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...