From 4a143d888ca1d4b3b46569b1784d04591ea1808d Mon Sep 17 00:00:00 2001 From: cpethrick <94074298+cpethrick@users.noreply.github.com> Date: Wed, 11 May 2022 14:24:23 -0400 Subject: [PATCH] Time Refinement Study using flow_solver (#138) * First commit for creating a time refinement study based on the flow_solver framework [does not build] * Flow solver runs 1D advection for one time-step size * Time refinement implemented and working * Writing convergence summary at each refinement * Bug fixes - test now runs as before * Adding exact solution [has bugs] * Bug fixes from previous commit * Results replicate 1D_explicit_ode_solver * Reorganized files and removed 1D_explicit_ode_solver test by removing directory from CMakeLists.txt * clean files/comments/naming * Renamed calculate_L2_error_at_final_time() --> calculate_L2_error_at_final_time_wrt_function() in anticipation of comparing to a reference solution' * Moved grid generation into straight_periodic_cube.cppw * Minor changes (comments, naming, minor functional changes) according to PR comments * Changed tolerance to 1E-13 rather than 0.5 dt for exiting the while loop in run_test() * Periodic1DFlow class now derives from PeriodicCubeFlow. Removed periodic_1D_flow.. * Update src/testing/time_refinement_study_advection.cpp Co-authored-by: Julien Brillon <43619715+jbrillon@users.noreply.github.com> * Renamed TimeRefinementStudyAdvection --> TimeRefinementStudy * Moved calculate_L2_error_at_final_time_wrt_function() into TimeRefinementStudy * Added error check that initial_time_step evenly divides final_time * Renaming files and small changes as per PR comments * Added time refinement study section to parameters; only passing FlowSolverParams in create_ExactSolutionFunction() * Fixed as per PR comments: createExactSolution now takes time argument; Periodic1DUnsteady renamed and moved into separate files; a few minor fixes * Bug fixes from merge * Update src/testing/time_refinement_study.cpp Co-authored-by: Julien Brillon <43619715+jbrillon@users.noreply.github.com> * Update src/mesh/grids/straight_periodic_cube.cpp Co-authored-by: Julien Brillon <43619715+jbrillon@users.noreply.github.com> * Changes as per PR comments: remove extra #include lines, initialize t in member initializer list, fix comments Co-authored-by: Julien Brillon <43619715+jbrillon@users.noreply.github.com> --- src/mesh/grids/straight_periodic_cube.cpp | 32 +++-- src/parameters/all_parameters.cpp | 7 +- src/parameters/all_parameters.h | 1 + src/parameters/parameters_flow_solver.cpp | 24 +++- src/parameters/parameters_flow_solver.h | 4 + src/physics/CMakeLists.txt | 1 + src/physics/exact_solutions/CMakeLists.txt | 22 +++ .../exact_solutions/exact_solution.cpp | 89 ++++++++++++ src/physics/exact_solutions/exact_solution.h | 83 +++++++++++ .../initial_conditions/initial_condition.cpp | 26 ++++ .../initial_conditions/initial_condition.h | 17 +++ src/testing/CMakeLists.txt | 7 +- src/testing/flow_solver.cpp | 8 +- src/testing/flow_solver.h | 3 +- .../flow_solver_case_base.cpp | 3 +- .../periodic_1D_unsteady.cpp | 50 +++++++ .../flow_solver_cases/periodic_1D_unsteady.h | 42 ++++++ .../flow_solver_cases/periodic_cube_flow.cpp | 5 +- .../flow_solver_cases/periodic_cube_flow.h | 2 +- src/testing/tests.cpp | 3 + src/testing/time_refinement_study.cpp | 132 ++++++++++++++++++ src/testing/time_refinement_study.h | 47 +++++++ .../CMakeLists.txt | 1 + .../flow_solver/CMakeLists.txt | 2 +- .../time_refinement_study/CMakeLists.txt | 15 ++ ...me_refinement_study_advection_explicit.prm | 44 ++++++ tests/unit_tests/CMakeLists.txt | 1 - 27 files changed, 648 insertions(+), 23 deletions(-) create mode 100644 src/physics/exact_solutions/CMakeLists.txt create mode 100644 src/physics/exact_solutions/exact_solution.cpp create mode 100644 src/physics/exact_solutions/exact_solution.h create mode 100644 src/testing/flow_solver_cases/periodic_1D_unsteady.cpp create mode 100644 src/testing/flow_solver_cases/periodic_1D_unsteady.h create mode 100644 src/testing/time_refinement_study.cpp create mode 100644 src/testing/time_refinement_study.h create mode 100644 tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt create mode 100644 tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm diff --git a/src/mesh/grids/straight_periodic_cube.cpp b/src/mesh/grids/straight_periodic_cube.cpp index 40d3482cb..7bbcd4fc5 100644 --- a/src/mesh/grids/straight_periodic_cube.cpp +++ b/src/mesh/grids/straight_periodic_cube.cpp @@ -24,24 +24,36 @@ void straight_periodic_cube(std::shared_ptr &grid, // Definition for each type of grid std::string grid_type_string; - if(dim==3) { + const bool colorize = true; + dealii::GridGenerator::hyper_cube(*grid, domain_left, domain_right, colorize); + if constexpr(dim == 1){ + grid_type_string = "Periodic 1D domain."; + std::vector > matched_pairs; + dealii::GridTools::collect_periodic_faces(*grid,0,1,0,matched_pairs); + grid->add_periodicity(matched_pairs); + }else if constexpr(dim==2) { + grid_type_string = "Doubly periodic square."; + std::vector::cell_iterator> > matched_pairs; + dealii::GridTools::collect_periodic_faces(*grid,0,1,0,matched_pairs); + dealii::GridTools::collect_periodic_faces(*grid,2,3,1,matched_pairs); + grid->add_periodicity(matched_pairs); + }else if constexpr(dim==3) { grid_type_string = "Triply periodic cube."; - const bool colorize = true; - dealii::GridGenerator::hyper_cube(*grid, domain_left, domain_right, colorize); std::vector::cell_iterator> > matched_pairs; dealii::GridTools::collect_periodic_faces(*grid,0,1,0,matched_pairs); dealii::GridTools::collect_periodic_faces(*grid,2,3,1,matched_pairs); dealii::GridTools::collect_periodic_faces(*grid,4,5,2,matched_pairs); grid->add_periodicity(matched_pairs); - grid->refine_global(number_of_refinements); - } - else { - std::cout << " ERROR: straight_periodic_cube() cannot be called for dim!=3 " << std::endl; - std::abort(); } + grid->refine_global(number_of_refinements); } -template void straight_periodic_cube<3, dealii::parallel::distributed::Triangulation<3>> (std::shared_ptr> &grid, const double domain_left, const double domain_right, const int number_of_cells_per_direction); +#if PHILIP_DIM==1 + template void straight_periodic_cube> (std::shared_ptr> &grid, const double domain_left, const double domain_right, const int number_of_cells_per_direction); +#endif +#if PHILIP_DIM!=1 + template void straight_periodic_cube> (std::shared_ptr> &grid, const double domain_left, const double domain_right, const int number_of_cells_per_direction); +#endif } // namespace Grids -} // namespace PHiLiP \ No newline at end of file +} // namespace PHiLiP diff --git a/src/parameters/all_parameters.cpp b/src/parameters/all_parameters.cpp index 2f57a2027..76bd7f13e 100644 --- a/src/parameters/all_parameters.cpp +++ b/src/parameters/all_parameters.cpp @@ -116,7 +116,8 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm) " flow_solver | " " dual_weighted_residual_mesh_adaptation | " " taylor_green_vortex_energy_check | " - " taylor_green_vortex_restart_check"), + " taylor_green_vortex_restart_check | " + " time_refinement_study"), "The type of test we want to solve. " "Choices are (only run control has been coded up for now)" " ."); + " taylor_green_vortex_restart_check | " + " time_refinement_study>."); prm.declare_entry("pde_type", "advection", dealii::Patterns::Selection( @@ -243,6 +245,7 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm) else if (test_string == "dual_weighted_residual_mesh_adaptation") { test_type = dual_weighted_residual_mesh_adaptation; } else if (test_string == "taylor_green_vortex_energy_check") { test_type = taylor_green_vortex_energy_check; } else if (test_string == "taylor_green_vortex_restart_check") { test_type = taylor_green_vortex_restart_check; } + else if (test_string == "time_refinement_study") { test_type = time_refinement_study; } const std::string pde_string = prm.get("pde_type"); if (pde_string == "advection") { diff --git a/src/parameters/all_parameters.h b/src/parameters/all_parameters.h index 27ab7b8b2..d22dbd9d3 100644 --- a/src/parameters/all_parameters.h +++ b/src/parameters/all_parameters.h @@ -131,6 +131,7 @@ class AllParameters dual_weighted_residual_mesh_adaptation, taylor_green_vortex_energy_check, taylor_green_vortex_restart_check, + time_refinement_study, }; TestType test_type; ///< Selected TestType from the input file. diff --git a/src/parameters/parameters_flow_solver.cpp b/src/parameters/parameters_flow_solver.cpp index 3d973a0bb..2b48ce76e 100644 --- a/src/parameters/parameters_flow_solver.cpp +++ b/src/parameters/parameters_flow_solver.cpp @@ -18,13 +18,15 @@ void FlowSolverParam::declare_parameters(dealii::ParameterHandler &prm) " taylor_green_vortex | " " burgers_viscous_snapshot | " " naca0012 | " - " burgers_rewienski_snapshot"), + " burgers_rewienski_snapshot | " + " advection_periodic"), "The type of flow we want to simulate. " "Choices are " " ."); + " burgers_rewienski_snapshot | " + " advection_periodic>."); prm.declare_entry("final_time", "1", dealii::Patterns::Double(0, dealii::Patterns::Double::max_double_value), @@ -94,6 +96,17 @@ void FlowSolverParam::declare_parameters(dealii::ParameterHandler &prm) " isothermal>."); } prm.leave_subsection(); + + prm.enter_subsection("time_refinement_study"); + { + prm.declare_entry("number_of_times_to_solve", "4", + dealii::Patterns::Integer(1, dealii::Patterns::Integer::max_int_value), + "Number of times to run the flow solver during a time refinement study."); + prm.declare_entry("refinement_ratio", "0.5", + dealii::Patterns::Double(0, 1.0), + "Ratio between a timestep size and the next in a time refinement study, 0 +#include "exact_solution.h" + +namespace PHiLiP { + +// ======================================================== +// ZERO -- Returns zero everywhere; used a placeholder when no exact solution is defined. +// ======================================================== +template +ExactSolutionFunction_Zero +::ExactSolutionFunction_Zero(double time_compare) + : ExactSolutionFunction() + , t(time_compare) +{ +} + +template +inline real ExactSolutionFunction_Zero +::value(const dealii::Point &/*point*/, const unsigned int /*istate*/) const +{ + real value = 0; + return value; +} + +// ======================================================== +// 1D SINE -- Exact solution for advection_explicit_time_study +// ======================================================== +template +ExactSolutionFunction_1DSine +::ExactSolutionFunction_1DSine (double time_compare) + : ExactSolutionFunction() + , t(time_compare) +{ +} + +template +inline real ExactSolutionFunction_1DSine +::value(const dealii::Point &point, const unsigned int /*istate*/) const +{ + double x_adv_speed = 1.0; + + real value = 0; + real pi = dealii::numbers::PI; + if(point[0] >= 0.0 && point[0] <= 2.0){ + value = sin(2*pi*(point[0] - x_adv_speed * t)/2.0); + } + return value; +} + + +//========================================================= +// FLOW SOLVER -- Exact Solution Base Class + Factory +//========================================================= +template +ExactSolutionFunction +::ExactSolutionFunction () + : dealii::Function(nstate) +{ + //do nothing +} + +template +std::shared_ptr> +ExactSolutionFactory::create_ExactSolutionFunction( + const Parameters::FlowSolverParam& flow_solver_parameters, + const double time_compare) +{ + // Get the flow case type + const FlowCaseEnum flow_type = flow_solver_parameters.flow_case_type; + if (flow_type == FlowCaseEnum::advection_periodic) { + if constexpr (dim==1 && nstate==dim) return std::make_shared > (time_compare); + } else { + // Select zero function if there is no exact solution defined + return std::make_shared> (time_compare); + } + return nullptr; +} + +template class ExactSolutionFunction ; +template class ExactSolutionFunction ; +template class ExactSolutionFactory ; +template class ExactSolutionFactory ; +template class ExactSolutionFunction_Zero ; +template class ExactSolutionFunction_Zero ; +template class ExactSolutionFunction_Zero ; +template class ExactSolutionFunction_Zero ; +template class ExactSolutionFunction_Zero ; + +} // PHiLiP namespace diff --git a/src/physics/exact_solutions/exact_solution.h b/src/physics/exact_solutions/exact_solution.h new file mode 100644 index 000000000..39f55ceaf --- /dev/null +++ b/src/physics/exact_solutions/exact_solution.h @@ -0,0 +1,83 @@ +#ifndef __EXACT_SOLUTION_H__ +#define __EXACT_SOLUTION_H__ + +// for the exact_solution function +#include +#include +#include "parameters/all_parameters.h" + +namespace PHiLiP { + +/// Initial condition function used to initialize a particular flow setup/case +template +class ExactSolutionFunction : public dealii::Function +{ +protected: + using dealii::Function::value; ///< dealii::Function we are templating on + +public: + /// Constructor + ExactSolutionFunction(); + /// Destructor + ~ExactSolutionFunction() {}; + + /// Value of the initial condition + virtual real value (const dealii::Point &point, const unsigned int istate = 0) const = 0; +}; + +/// Initial Condition Function: Zero Function; used as a placeholder when there is no exact solution +template +class ExactSolutionFunction_Zero + : public ExactSolutionFunction +{ +protected: + using dealii::Function::value; ///< dealii::Function we are templating on + +public: + /// Constructor for InitialConditionFunction_Zero + ExactSolutionFunction_Zero (double time_compare); + + /// Time at which to compute the exact solution + double t; + + /// Value of initial condition + real value (const dealii::Point &point, const unsigned int istate = 0) const override; +}; + +/// Initial Condition Function: 1D Sine Function; used for temporal convergence +template +class ExactSolutionFunction_1DSine + : public ExactSolutionFunction +{ +protected: + using dealii::Function::value; ///< dealii::Function we are templating on + +public: + /// Constructor for InitialConditionFunction_1DSine + /** Calls the Function(const unsigned int n_components) constructor in deal.II*/ + ExactSolutionFunction_1DSine (double time_compare); + + /// Time at which to compute the exact solution + double t; + + /// Value of initial condition + real value (const dealii::Point &point, const unsigned int istate = 0) const override; +}; + + +/// Initial condition function factory +template +class ExactSolutionFactory +{ +protected: + /// Enumeration of all flow solver initial conditions types defined in the Parameters class + using FlowCaseEnum = Parameters::FlowSolverParam::FlowCaseType; + +public: + /// Construct InitialConditionFunction object from global parameter file + static std::shared_ptr> + create_ExactSolutionFunction(const Parameters::FlowSolverParam& flow_solver_parameters, const double time_compare); +}; + +} // PHiLiP namespace +#endif diff --git a/src/physics/initial_conditions/initial_condition.cpp b/src/physics/initial_conditions/initial_condition.cpp index 3e0a81202..faea30dab 100644 --- a/src/physics/initial_conditions/initial_condition.cpp +++ b/src/physics/initial_conditions/initial_condition.cpp @@ -158,6 +158,30 @@ ::value(const dealii::Point &point, const unsigned int /*istate*/) con } +// ======================================================== +// 1D SINE -- Initial Condition for advection_explicit_time_study +// ======================================================== +template +InitialConditionFunction_1DSine +::InitialConditionFunction_1DSine () + : InitialConditionFunction() +{ + // Nothing to do here yet +} + +template +inline real InitialConditionFunction_1DSine +::value(const dealii::Point &point, const unsigned int /*istate*/) const +{ + real value = 0; + real pi = dealii::numbers::PI; + if(point[0] >= 0.0 && point[0] <= 2.0){ + value = sin(2*pi*point[0]/2.0); + } + return value; +} + + //========================================================= // FLOW SOLVER -- Initial Condition Base Class + Factory //========================================================= @@ -204,6 +228,8 @@ InitialConditionFactory::create_InitialConditionFunction( param->euler_param.side_slip_angle); return std::make_shared>(euler_physics_double); } + } else if (flow_type == FlowCaseEnum::advection_periodic) { + if constexpr (dim==1 && nstate==dim) return std::make_shared > (); } else { std::cout << "Invalid Flow Case Type. You probably forgot to add it to the list of flow cases in initial_condition.cpp" << std::endl; std::abort(); diff --git a/src/physics/initial_conditions/initial_condition.h b/src/physics/initial_conditions/initial_condition.h index 788de099f..3dbb285ab 100644 --- a/src/physics/initial_conditions/initial_condition.h +++ b/src/physics/initial_conditions/initial_condition.h @@ -151,6 +151,23 @@ class InitialConditionFunction_BurgersViscous: public InitialConditionFunction &point, const unsigned int istate = 0) const override; }; +/// Initial Condition Function: 1D Sine Function; used for temporal convergence +template +class InitialConditionFunction_1DSine + : public InitialConditionFunction +{ +protected: + using dealii::Function::value; ///< dealii::Function we are templating on + +public: + /// Constructor for InitialConditionFunction_1DSine + InitialConditionFunction_1DSine (); + + /// Value of initial condition + real value (const dealii::Point &point, const unsigned int istate = 0) const override; +}; + + /// Initial condition function factory template class InitialConditionFactory diff --git a/src/testing/CMakeLists.txt b/src/testing/CMakeLists.txt index 9d6c47cef..57d582eb9 100644 --- a/src/testing/CMakeLists.txt +++ b/src/testing/CMakeLists.txt @@ -25,13 +25,15 @@ set(TESTS_SOURCE flow_solver_cases/flow_solver_case_base.cpp flow_solver_cases/periodic_cube_flow.cpp flow_solver_cases/periodic_turbulence.cpp + flow_solver_cases/periodic_1D_unsteady.cpp flow_solver_cases/1D_burgers_rewienski_snapshot.cpp flow_solver_cases/1d_burgers_viscous_snapshot.cpp flow_solver.cpp dual_weighted_residual_mesh_adaptation.cpp taylor_green_vortex_energy_check.cpp taylor_green_vortex_restart_check.cpp - flow_solver_cases/naca0012.cpp) + flow_solver_cases/naca0012.cpp + time_refinement_study.cpp) foreach(dim RANGE 1 3) # Output library @@ -52,6 +54,7 @@ foreach(dim RANGE 1 3) string(CONCAT GridRefinementLib GridRefinement_${dim}D) string(CONCAT PODLib POD_${dim}D) string(CONCAT InitialConditionsLib InitialConditions_${dim}D) + string(CONCAT ExactSolutionsLib ExactSolutions_${dim}D) target_link_libraries(${TestsLib} ${GridsLib}) target_link_libraries(${TestsLib} ${NumericalFluxLib}) target_link_libraries(${TestsLib} ${PhysicsLib}) @@ -63,6 +66,7 @@ foreach(dim RANGE 1 3) target_link_libraries(${TestsLib} ${OptimizationLib}) target_link_libraries(${TestsLib} ${PODLib}) target_link_libraries(${TestsLib} ${InitialConditionsLib}) + target_link_libraries(${TestsLib} ${ExactSolutionsLib}) # Setup target with deal.II if(NOT DOC_ONLY) DEAL_II_SETUP_TARGET(${TestsLib}) @@ -80,4 +84,5 @@ foreach(dim RANGE 1 3) unset(OptimizationLib) unset(PODLib) unset(InitialConditionsLib) + unset(ExactSolutionsLib) endforeach() diff --git a/src/testing/flow_solver.cpp b/src/testing/flow_solver.cpp index 74b410de3..534b99d28 100644 --- a/src/testing/flow_solver.cpp +++ b/src/testing/flow_solver.cpp @@ -7,6 +7,7 @@ #include #include + namespace PHiLiP { namespace Tests { @@ -371,7 +372,7 @@ int FlowSolver::run_test() const // Time advancement loop with on-the-fly post-processing //---------------------------------------------------- pcout << "Advancing solution in time... " << std::endl; - while(ode_solver->current_time < final_time) + while((ode_solver->current_time) < (final_time - 1E-13)) //comparing to 1E-13 to avoid taking an extra timestep { // advance solution ode_solver->step_in_time(constant_time_step,false); // pseudotime==false @@ -462,6 +463,11 @@ ::create_FlowSolver(const Parameters::AllParameters *const parameters_input, std::shared_ptr> flow_solver_case = std::make_shared>(parameters_input); return std::make_unique>(parameters_input, flow_solver_case, parameter_handler_input); } + } else if (flow_type == FlowCaseEnum::advection_periodic){ + if constexpr (dim==1 && nstate==dim){ + std::shared_ptr> flow_solver_case = std::make_shared>(parameters_input); + return std::make_unique>(parameters_input, flow_solver_case, parameter_handler_input); + } } else { std::cout << "Invalid flow case. You probably forgot to add it to the list of flow cases in flow_solver.cpp" << std::endl; std::abort(); diff --git a/src/testing/flow_solver.h b/src/testing/flow_solver.h index 95537dd62..1b7835a8f 100644 --- a/src/testing/flow_solver.h +++ b/src/testing/flow_solver.h @@ -26,6 +26,7 @@ #include "ode_solver/ode_solver_factory.h" #include "flow_solver_cases/periodic_cube_flow.h" #include "flow_solver_cases/periodic_turbulence.h" +#include "flow_solver_cases/periodic_1D_unsteady.h" #include "flow_solver_cases/1D_burgers_rewienski_snapshot.h" #include "flow_solver_cases/1d_burgers_viscous_snapshot.h" #include "flow_solver_cases/naca0012.h" @@ -127,4 +128,4 @@ class FlowSolverFactory } // Tests namespace } // PHiLiP namespace -#endif \ No newline at end of file +#endif diff --git a/src/testing/flow_solver_cases/flow_solver_case_base.cpp b/src/testing/flow_solver_cases/flow_solver_case_base.cpp index 0043ad9a3..eec6e5956 100644 --- a/src/testing/flow_solver_cases/flow_solver_case_base.cpp +++ b/src/testing/flow_solver_cases/flow_solver_case_base.cpp @@ -46,6 +46,7 @@ std::string FlowSolverCaseBase::get_flow_case_string() const if (flow_case_type == FlowCaseEnum::burgers_viscous_snapshot) {flow_case_string = "burgers_viscous_snapshot";} if (flow_case_type == FlowCaseEnum::burgers_rewienski_snapshot) {flow_case_string = "burgers_rewienski_snapshot";} if (flow_case_type == FlowCaseEnum::naca0012) {flow_case_string = "naca0012";} + if (flow_case_type == FlowCaseEnum::advection_periodic) {flow_case_string = "advection_periodic";} return flow_case_string; } @@ -117,4 +118,4 @@ template class FlowSolverCaseBase; #endif } -} \ No newline at end of file +} diff --git a/src/testing/flow_solver_cases/periodic_1D_unsteady.cpp b/src/testing/flow_solver_cases/periodic_1D_unsteady.cpp new file mode 100644 index 000000000..7434430e7 --- /dev/null +++ b/src/testing/flow_solver_cases/periodic_1D_unsteady.cpp @@ -0,0 +1,50 @@ +#include "periodic_1D_unsteady.h" +#include +#include +#include +#include +#include +#include +#include +#include "physics/physics_factory.h" +#include "dg/dg.h" +#include "mesh/grids/straight_periodic_cube.hpp" + +namespace PHiLiP { + +namespace Tests { + +//========================================================= +// PERIODIC 1D DOMAIN FOR UNSTEADY CALCULATIONS +//========================================================= + +template +Periodic1DUnsteady::Periodic1DUnsteady(const PHiLiP::Parameters::AllParameters *const parameters_input) + : PeriodicCubeFlow(parameters_input) +{ +} + +template +void Periodic1DUnsteady::compute_unsteady_data_and_write_to_table( + const unsigned int current_iteration, + const double current_time, + const std::shared_ptr >/* dg */, + const std::shared_ptr /* unsteady_data_table */) +{ + double dt = this->all_param.ode_solver_param.initial_time_step; + int output_solution_every_n_iterations = (int) (this->all_param.ode_solver_param.output_solution_every_dt_time_intervals/dt); + + if ((current_iteration % output_solution_every_n_iterations) == 0){ + this->pcout << " Iter: " << current_iteration + << " Time: " << current_time + << std::endl; + } +} + +#if PHILIP_DIM==1 +template class Periodic1DUnsteady ; +#endif + +} // Tests namespace +} // PHiLiP namespace + diff --git a/src/testing/flow_solver_cases/periodic_1D_unsteady.h b/src/testing/flow_solver_cases/periodic_1D_unsteady.h new file mode 100644 index 000000000..b53a258bc --- /dev/null +++ b/src/testing/flow_solver_cases/periodic_1D_unsteady.h @@ -0,0 +1,42 @@ +#ifndef __PERIODIC_1D_UNSTEADY_H__ +#define __PERIODIC_1D_UNSTEADY_H__ + +// for FlowSolver class: +#include "periodic_cube_flow.h" +#include "dg/dg.h" +#include "physics/physics.h" +#include "parameters/all_parameters.h" +#include + +namespace PHiLiP { +namespace Tests { + +template +class Periodic1DUnsteady : public PeriodicCubeFlow +{ +public: + + /// Constructor. + Periodic1DUnsteady(const Parameters::AllParameters *const parameters_input); + + /// Destructor + ~Periodic1DUnsteady() {}; +protected: + + /// Compute the desired unsteady data and write it to a table + /** Currently only prints to console. + * In the future, will be modified to also write to table if there is unsteady + * data of interest + */ + void compute_unsteady_data_and_write_to_table( + const unsigned int current_iteration, + const double current_time, + const std::shared_ptr > dg, + const std::shared_ptr unsteady_data_table) override; + +}; + + +} // Tests namespace +} // PHiLiP namespace +#endif diff --git a/src/testing/flow_solver_cases/periodic_cube_flow.cpp b/src/testing/flow_solver_cases/periodic_cube_flow.cpp index 44ea35355..222c7618b 100644 --- a/src/testing/flow_solver_cases/periodic_cube_flow.cpp +++ b/src/testing/flow_solver_cases/periodic_cube_flow.cpp @@ -27,7 +27,7 @@ std::shared_ptr PeriodicCubeFlow::generate_grid() con this->mpi_communicator #endif ); - Grids::straight_periodic_cube>(grid, domain_left, domain_right, number_of_cells_per_direction); + Grids::straight_periodic_cube(grid, domain_left, domain_right, number_of_cells_per_direction); return grid; } @@ -55,8 +55,9 @@ void PeriodicCubeFlow::display_additional_flow_case_specific_paramet #if PHILIP_DIM==3 template class PeriodicCubeFlow ; +#elif PHILIP_DIM==1 +template class PeriodicCubeFlow ; #endif - } // Tests namespace } // PHiLiP namespace diff --git a/src/testing/flow_solver_cases/periodic_cube_flow.h b/src/testing/flow_solver_cases/periodic_cube_flow.h index 3ab0fdbcf..14e995654 100644 --- a/src/testing/flow_solver_cases/periodic_cube_flow.h +++ b/src/testing/flow_solver_cases/periodic_cube_flow.h @@ -40,4 +40,4 @@ class PeriodicCubeFlow : public FlowSolverCaseBase } // Tests namespace } // PHiLiP namespace -#endif \ No newline at end of file +#endif diff --git a/src/testing/tests.cpp b/src/testing/tests.cpp index 11d264911..654626da4 100644 --- a/src/testing/tests.cpp +++ b/src/testing/tests.cpp @@ -34,6 +34,7 @@ #include "dual_weighted_residual_mesh_adaptation.h" #include "taylor_green_vortex_energy_check.h" #include "taylor_green_vortex_restart_check.h" +#include "time_refinement_study.h" namespace PHiLiP { namespace Tests { @@ -169,6 +170,8 @@ ::select_test(const AllParam *const parameters_input, if constexpr (dim==3 && nstate==dim+2) return std::make_unique>(parameters_input,parameter_handler_input); } else if(test_type == Test_enum::taylor_green_vortex_restart_check) { if constexpr (dim==3 && nstate==dim+2) return std::make_unique>(parameters_input,parameter_handler_input); + } else if(test_type == Test_enum::time_refinement_study) { + if constexpr (dim==1 && nstate==1) return std::make_unique>(parameters_input, parameter_handler_input); } else { std::cout << "Invalid test. You probably forgot to add it to the list of tests in tests.cpp" << std::endl; std::abort(); diff --git a/src/testing/time_refinement_study.cpp b/src/testing/time_refinement_study.cpp new file mode 100644 index 000000000..7abd968f8 --- /dev/null +++ b/src/testing/time_refinement_study.cpp @@ -0,0 +1,132 @@ +#include "time_refinement_study.h" +#include "flow_solver.h" +#include "flow_solver_cases/periodic_1D_unsteady.h" +#include "physics/exact_solutions/exact_solution.h" +#include "cmath" + +namespace PHiLiP { +namespace Tests { + +template +TimeRefinementStudy::TimeRefinementStudy( + const PHiLiP::Parameters::AllParameters *const parameters_input, + const dealii::ParameterHandler ¶meter_handler_input) + : TestsBase::TestsBase(parameters_input), + parameter_handler(parameter_handler_input), + n_time_calculations(parameters_input->flow_solver_param.number_of_times_to_solve), + refine_ratio(parameters_input->flow_solver_param.refinement_ratio) +{} + + + +template +Parameters::AllParameters TimeRefinementStudy::reinit_params_and_refine_timestep(int refinement) const +{ + PHiLiP::Parameters::AllParameters parameters = *(this->all_parameters); + + parameters.ode_solver_param.initial_time_step *= pow(refine_ratio,refinement); + + return parameters; +} + +template +double TimeRefinementStudy::calculate_L2_error_at_final_time_wrt_function(std::shared_ptr> dg, const Parameters::AllParameters parameters, double final_time) const +{ + + //generate exact solution at final time + std::shared_ptr> exact_solution_function; + exact_solution_function = ExactSolutionFactory::create_ExactSolutionFunction(parameters.flow_solver_param, final_time); + int poly_degree = parameters.grid_refinement_study_param.poly_degree; + dealii::Vector difference_per_cell(dg->solution.size()); + + dealii::VectorTools::integrate_difference(dg->dof_handler, + dg->solution, + *exact_solution_function, + difference_per_cell, + dealii::QGauss(poly_degree+10), //overintegrating by 10 + dealii::VectorTools::L2_norm); + + double L2_error = dealii::VectorTools::compute_global_error(*dg->triangulation, + difference_per_cell, + dealii::VectorTools::L2_norm); + return L2_error; + +} + + +template +int TimeRefinementStudy::run_test() const +{ + + double final_time = this->all_parameters->flow_solver_param.final_time; + double initial_time_step = this->all_parameters->ode_solver_param.initial_time_step; + int n_steps = round(final_time/initial_time_step); + if (n_steps * initial_time_step != final_time){ + pcout << "Error: final_time is not evenly divisible by initial_time_step!" << std::endl + << "Remainder is " << fmod(final_time, initial_time_step) + << ". Modify parameters to run this test." << std::endl; + std::abort(); + } + + int testfail = 0; + double expected_order =(double) this->all_parameters->ode_solver_param.runge_kutta_order; + double order_tolerance = 0.1; + + dealii::ConvergenceTable convergence_table; + double L2_error_old = 0; + double L2_error_conv_rate=0; + + for (int refinement = 0; refinement < n_time_calculations; ++refinement){ + + pcout << "\n\n---------------------------------------------" << std::endl; + pcout << "Refinement number " << refinement << " of " << n_time_calculations - 1 << std::endl; + pcout << "---------------------------------------------" << std::endl; + + const Parameters::AllParameters params = reinit_params_and_refine_timestep(refinement); + std::unique_ptr> flow_solver = FlowSolverFactory::create_FlowSolver(¶ms, parameter_handler); + static_cast(flow_solver->run_test()); + + //check L2 error + double L2_error = calculate_L2_error_at_final_time_wrt_function(flow_solver->dg, params, flow_solver->ode_solver->current_time); + pcout << "Computed error is " << L2_error << std::endl; + + const double dt = params.ode_solver_param.initial_time_step; + convergence_table.add_value("refinement", refinement); + convergence_table.add_value("dt", dt ); + convergence_table.set_precision("dt", 16); + convergence_table.set_scientific("dt", true); + convergence_table.add_value("L2_error",L2_error); + convergence_table.set_precision("L2_error", 16); + convergence_table.evaluate_convergence_rates("L2_error", "dt", dealii::ConvergenceTable::reduction_rate_log2, 1); + + //Checking convergence order + if (refinement > 0) { + L2_error_conv_rate = -log(L2_error_old/L2_error)/log(refine_ratio); + pcout << "Order at " << refinement << " is " << L2_error_conv_rate << std::endl; + if (abs(L2_error_conv_rate - expected_order) > order_tolerance){ + testfail = 1; + pcout << "Expected convergence order was not reached at refinement " << refinement <; +#endif +} // Tests namespace +} // PHiLiP namespace diff --git a/src/testing/time_refinement_study.h b/src/testing/time_refinement_study.h new file mode 100644 index 000000000..109def64e --- /dev/null +++ b/src/testing/time_refinement_study.h @@ -0,0 +1,47 @@ +#ifndef __TIME_REFINEMENT_STUDY_ADVECTION__ +#define __TIME_REFINEMENT_STUDY_ADVECTION__ + +#include "tests.h" +#include +#include "dg/dg.h" + +namespace PHiLiP { +namespace Tests { + +/// Advection time refinement study +template +class TimeRefinementStudy: public TestsBase +{ +public: + /// Constructor + TimeRefinementStudy( + const Parameters::AllParameters *const parameters_input, + const dealii::ParameterHandler ¶meter_handler_input); + + /// Destructor + ~TimeRefinementStudy() {}; + + /// Parameter handler for storing the .prm file being ran + const dealii::ParameterHandler ¶meter_handler; + + /// Run test + int run_test () const override; +protected: + /// Number of times to solve for convergence summary + const int n_time_calculations; + + /// Ratio to refine by + const double refine_ratio; + + /// Calculate L2 error at the final time in the passed parameters + double calculate_L2_error_at_final_time_wrt_function(std::shared_ptr> dg,const Parameters::AllParameters parameters, double final_time) const; + + /// Reinitialize parameters while refining the timestep. Necessary because all_parameters is constant. + Parameters::AllParameters reinit_params_and_refine_timestep(int refinement) const; + +}; + +} // End of Tests namespace +} // End of PHiLiP namespace + +#endif diff --git a/tests/integration_tests_control_files/CMakeLists.txt b/tests/integration_tests_control_files/CMakeLists.txt index bed73ac63..85d5623ec 100644 --- a/tests/integration_tests_control_files/CMakeLists.txt +++ b/tests/integration_tests_control_files/CMakeLists.txt @@ -17,3 +17,4 @@ add_subdirectory(1d_shock) add_subdirectory(reduced_order) add_subdirectory(flow_solver) add_subdirectory(taylor_green_vortex_integration) +add_subdirectory(time_refinement_study) diff --git a/tests/integration_tests_control_files/flow_solver/CMakeLists.txt b/tests/integration_tests_control_files/flow_solver/CMakeLists.txt index 4845c3d5a..528fd50fe 100644 --- a/tests/integration_tests_control_files/flow_solver/CMakeLists.txt +++ b/tests/integration_tests_control_files/flow_solver/CMakeLists.txt @@ -33,4 +33,4 @@ set(TEST_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}) # COMMAND mpirun -np ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/viscous_taylor_green_vortex.prm # WORKING_DIRECTORY ${TEST_OUTPUT_DIR} #) -# ---------------------------------------- \ No newline at end of file +# ---------------------------------------- diff --git a/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt b/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt new file mode 100644 index 000000000..fd012ebd7 --- /dev/null +++ b/tests/integration_tests_control_files/time_refinement_study/CMakeLists.txt @@ -0,0 +1,15 @@ +set(TEST_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}) + +# ======================================= +# Time Study (Linear Advection Explicit RK) +# ======================================= +# ---------------------------------------- +# - details: [if needed] +# ---------------------------------------- +configure_file(time_refinement_study_advection_explicit.prm time_refinement_study_advection_explicit.prm COPYONLY) +add_test( + NAME 1D_TIME_REFINEMENT_STUDY_ADVECTION_EXPLICIT + COMMAND mpirun -np 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/time_refinement_study_advection_explicit.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) +# ---------------------------------------- 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 new file mode 100644 index 000000000..b75343a32 --- /dev/null +++ b/tests/integration_tests_control_files/time_refinement_study/time_refinement_study_advection_explicit.prm @@ -0,0 +1,44 @@ +# 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 = explicit + set output_solution_every_dt_time_intervals = 0.1 + set initial_time_step = 2.5E-3 + subsection explicit solver options + set runge_kutta_order = 3 + end +end + +subsection manufactured solution convergence study + # advection speed + set advection_0 = 1.0 + set advection_1 = 0.0 +end + +# polynomial order and number of cells per direction (i.e. grid_size) +subsection grid refinement study + set poly_degree = 5 + #set num_refinements = 5 #using grid_size instead + set grid_size = 32 + set grid_left = 0.0 + set grid_right = 2.0 +end + +subsection flow_solver + set flow_case_type = advection_periodic + set final_time = 1.0 + subsection time_refinement_study + set number_of_times_to_solve = 4 + set refinement_ratio = 0.5 + end +end diff --git a/tests/unit_tests/CMakeLists.txt b/tests/unit_tests/CMakeLists.txt index 1c5991f2b..77ba7544b 100644 --- a/tests/unit_tests/CMakeLists.txt +++ b/tests/unit_tests/CMakeLists.txt @@ -9,4 +9,3 @@ add_subdirectory(msh_output) add_subdirectory(navier_stokes_unit_test) add_subdirectory(optimization) add_subdirectory(operator_tests) -add_subdirectory(ode_solver_unit_test)