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.<cpp/h>.

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

Author: cpethrick

src/mesh/grids/straight_periodic_cube.cpp

@@ -24,24 +24,36 @@ void straight_periodic_cube(std::shared_ptr<TriangulationType> &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<dealii::GridTools::PeriodicFacePair<typename TriangulationType::cell_iterator> > 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<dealii::GridTools::PeriodicFacePair<typename dealii::Triangulation<dim>::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<dealii::GridTools::PeriodicFacePair<typename dealii::Triangulation<dim>::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<dealii::parallel::distributed::Triangulation<3>> &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<PHILIP_DIM, dealii::Triangulation<PHILIP_DIM>> (std::shared_ptr<dealii::Triangulation<PHILIP_DIM>> &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<PHILIP_DIM, dealii::parallel::distributed::Triangulation<PHILIP_DIM>> (std::shared_ptr<dealii::parallel::distributed::Triangulation<PHILIP_DIM>> &grid, const double domain_left, const double domain_right, const int number_of_cells_per_direction);
+#endif
 
 } // namespace Grids
-} // namespace PHiLiP
+} // namespace PHiLiP

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)" 
                       " <run_control | " 
@@ -143,7 +144,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>.");
 
     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") {

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.
 

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 "
                           " <taylor_green_vortex | "
                           " burgers_viscous_snapshot | "
                           " naca0012 | "
-                          " burgers_rewienski_snapshot>.");
+                          " 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<r<1.");
+        }
+        prm.leave_subsection();
     }
     prm.leave_subsection();
 }
@@ -107,6 +120,7 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm)
         else if (flow_case_type_string == "burgers_viscous_snapshot")   {flow_case_type = burgers_viscous_snapshot;}
         else if (flow_case_type_string == "burgers_rewienski_snapshot") {flow_case_type = burgers_rewienski_snapshot;}
         else if (flow_case_type_string == "naca0012")                   {flow_case_type = naca0012;}
+        else if (flow_case_type_string == "advection_periodic")         {flow_case_type = advection_periodic;}
 
         final_time = prm.get_double("final_time");
         courant_friedrich_lewy_number = prm.get_double("courant_friedrich_lewy_number");
@@ -131,6 +145,12 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm)
             else if (density_initial_condition_type_string == "isothermal") {density_initial_condition_type = isothermal;}
         }
         prm.leave_subsection();
+        prm.enter_subsection("time_refinement_study");
+        {
+            number_of_times_to_solve = prm.get_integer("number_of_times_to_solve");
+            refinement_ratio = prm.get_double("refinement_ratio");
+        }
+        prm.leave_subsection();
     }
     prm.leave_subsection();
 }

src/parameters/parameters_flow_solver.h

@@ -20,6 +20,7 @@ class FlowSolverParam
         burgers_viscous_snapshot,
         naca0012,
         burgers_rewienski_snapshot,
+        advection_periodic
         };
     FlowCaseType flow_case_type; ///< Selected FlowCaseType from the input file
 
@@ -62,6 +63,9 @@ class FlowSolverParam
         };
     /// Selected DensityInitialConditionType from the input file
     DensityInitialConditionType density_initial_condition_type;
+  
+    int number_of_times_to_solve; ///<For time refinement study, number of times to run the calculation
+    double refinement_ratio; ///<For time refinement study, ratio of next timestep size to current one, 0<r<1
 
     /// Declares the possible variables and sets the defaults.
     static void declare_parameters (dealii::ParameterHandler &prm);

src/physics/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_subdirectory(initial_conditions)
+add_subdirectory(exact_solutions)
 
 SET(SOURCE
     physics_factory.cpp

src/physics/exact_solutions/CMakeLists.txt

@@ -0,0 +1,22 @@
+set(SOURCE
+    exact_solution.cpp
+    )
+
+foreach(dim RANGE 1 3)
+    # Output library
+    string(CONCAT ExactSolutionsLib ExactSolutions_${dim}D)
+    add_library(${ExactSolutionsLib} STATIC ${SOURCE})
+    # Link with PhysicsLib
+    string(CONCAT PhysicsLib Physics_${dim}D)
+    target_link_libraries(${ExactSolutionsLib} ${PhysicsLib})
+
+    target_compile_definitions(${ExactSolutionsLib} PRIVATE PHILIP_DIM=${dim})
+    unset(PhysicsLib)
+
+    # Setup target with deal.II
+    if(NOT DOC_ONLY)
+        DEAL_II_SETUP_TARGET(${ExactSolutionsLib})
+    endif()
+
+    unset(ExactSolutionsLib)
+endforeach()

src/physics/exact_solutions/exact_solution.cpp

@@ -0,0 +1,89 @@
+#include <deal.II/base/function.h>
+#include "exact_solution.h"
+
+namespace PHiLiP {
+
+// ========================================================
+// ZERO -- Returns zero everywhere; used a placeholder when no exact solution is defined.
+// ========================================================
+template <int dim, int nstate, typename real>
+ExactSolutionFunction_Zero<dim,nstate,real>
+::ExactSolutionFunction_Zero(double time_compare)
+        : ExactSolutionFunction<dim,nstate,real>()
+        , t(time_compare)
+{
+}
+
+template <int dim, int nstate, typename real>
+inline real ExactSolutionFunction_Zero<dim,nstate,real>
+::value(const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/) const
+{
+    real value = 0;
+    return value;
+}
+
+// ========================================================
+// 1D SINE -- Exact solution for advection_explicit_time_study
+// ========================================================
+template <int dim, int nstate, typename real>
+ExactSolutionFunction_1DSine<dim,nstate,real>
+::ExactSolutionFunction_1DSine (double time_compare)
+        : ExactSolutionFunction<dim,nstate,real>()
+        , t(time_compare)
+{
+}
+
+template <int dim, int nstate, typename real>
+inline real ExactSolutionFunction_1DSine<dim,nstate,real>
+::value(const dealii::Point<dim,real> &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 <int dim, int nstate, typename real>
+ExactSolutionFunction<dim,nstate,real>
+::ExactSolutionFunction ()
+    : dealii::Function<dim,real>(nstate)
+{
+    //do nothing
+}
+
+template <int dim, int nstate, typename real>
+std::shared_ptr<ExactSolutionFunction<dim, nstate, real>>
+ExactSolutionFactory<dim,nstate, real>::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<ExactSolutionFunction_1DSine<dim,nstate,real> > (time_compare);
+    } else {
+        // Select zero function if there is no exact solution defined
+        return std::make_shared<ExactSolutionFunction_Zero<dim,nstate,real>> (time_compare);
+    }
+    return nullptr;
+}
+
+template class ExactSolutionFunction <PHILIP_DIM,PHILIP_DIM, double>;
+template class ExactSolutionFunction <PHILIP_DIM,PHILIP_DIM+2, double>;
+template class ExactSolutionFactory <PHILIP_DIM, PHILIP_DIM+2, double>;
+template class ExactSolutionFactory <PHILIP_DIM, PHILIP_DIM, double>;
+template class ExactSolutionFunction_Zero <PHILIP_DIM,1, double>;
+template class ExactSolutionFunction_Zero <PHILIP_DIM,2, double>;
+template class ExactSolutionFunction_Zero <PHILIP_DIM,3, double>;
+template class ExactSolutionFunction_Zero <PHILIP_DIM,4, double>;
+template class ExactSolutionFunction_Zero <PHILIP_DIM,5, double>;
+
+} // PHiLiP namespace

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 <deal.II/lac/vector.h>
+#include <deal.II/base/function.h>
+#include "parameters/all_parameters.h"
+
+namespace PHiLiP {
+
+/// Initial condition function used to initialize a particular flow setup/case
+template <int dim, int nstate, typename real>
+class ExactSolutionFunction : public dealii::Function<dim,real>
+{
+protected:
+    using dealii::Function<dim,real>::value; ///< dealii::Function we are templating on
+
+public:
+    /// Constructor
+    ExactSolutionFunction();
+    /// Destructor
+    ~ExactSolutionFunction() {};
+
+    /// Value of the initial condition
+    virtual real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const = 0;
+};
+
+/// Initial Condition Function: Zero Function; used as a placeholder when there is no exact solution
+template <int dim, int nstate, typename real>
+class ExactSolutionFunction_Zero
+        : public ExactSolutionFunction<dim,nstate,real>
+{
+protected:
+    using dealii::Function<dim,real>::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<dim,real> &point, const unsigned int istate = 0) const override;
+};
+
+/// Initial Condition Function: 1D Sine Function; used for temporal convergence
+template <int dim, int nstate, typename real>
+class ExactSolutionFunction_1DSine
+        : public ExactSolutionFunction<dim,nstate,real>
+{
+protected:
+    using dealii::Function<dim,real>::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<dim,real> &point, const unsigned int istate = 0) const override;
+};
+
+
+/// Initial condition function factory
+template <int dim, int nstate, typename real>
+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<ExactSolutionFunction<dim,nstate,real>>
+        create_ExactSolutionFunction(const Parameters::FlowSolverParam& flow_solver_parameters, const double time_compare);
+};
+
+} // PHiLiP namespace
+#endif