Skip to content

Commit

Permalink
Time Refinement Study using flow_solver (#138)
Browse files Browse the repository at this point in the history
* 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>
  • Loading branch information
cpethrick and jbrillon authored May 11, 2022
1 parent b13f5ff commit 4a143d8
Show file tree
Hide file tree
Showing 27 changed files with 648 additions and 23 deletions.
32 changes: 22 additions & 10 deletions src/mesh/grids/straight_periodic_cube.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
7 changes: 5 additions & 2 deletions src/parameters/all_parameters.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 | "
Expand All @@ -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(
Expand Down Expand Up @@ -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") {
Expand Down
1 change: 1 addition & 0 deletions src/parameters/all_parameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
24 changes: 22 additions & 2 deletions src/parameters/parameters_flow_solver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down Expand Up @@ -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();
}
Expand All @@ -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");
Expand All @@ -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();
}
Expand Down
4 changes: 4 additions & 0 deletions src/parameters/parameters_flow_solver.h
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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);
Expand Down
1 change: 1 addition & 0 deletions src/physics/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
add_subdirectory(initial_conditions)
add_subdirectory(exact_solutions)

SET(SOURCE
physics_factory.cpp
Expand Down
22 changes: 22 additions & 0 deletions src/physics/exact_solutions/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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()
89 changes: 89 additions & 0 deletions src/physics/exact_solutions/exact_solution.cpp
Original file line number Diff line number Diff line change
@@ -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
83 changes: 83 additions & 0 deletions src/physics/exact_solutions/exact_solution.h
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 4a143d8

Please sign in to comment.