From fd24c81207ef418649eb7d1b719c54c86c7554c2 Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Tue, 17 Sep 2019 12:13:53 -0400 Subject: [PATCH 01/10] Update README.md --- README.md | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 3ab8345d6..409ac0a99 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,10 @@ ## Code Description - Code uses deal.II library as the backbone (https://www.dealii.org/) -- Supported Partial Differential Equations: Convection-diffusion, Euler, TODO: Navier-Stokes. -- Supported convective numerical fluxes: Lax-Friedrichs, Roe (Harten's entropy fix) for Euler +- Math supporting this code can be viewed in this **very rough draft in progress** [Overleaf document](https://www.overleaf.com/read/mytvbbbbyqnj). +- Supports weak and strong (InProgress) form of discontinuous Galerkin (DG), and flux reconstruction (FR) (InProgress) +- Supported Partial Differential Equations: Linear advection, diffusion, convection-diffusion, Burgers, Euler, TODO: Navier-Stokes. +- Supported convective numerical fluxes: Lax-Friedrichs, Roe (Harten's entropy fix) for Euler, InProgress: Split-Form - Supported diffusive numerical fluxes: Symmetric Interior Penalty - Supported elements: LINEs, QUADs, HEXs - Supported refinements: h (size) or p (order). @@ -95,11 +97,11 @@ GDB$ quit ### Parallel debugging -If the error only occurs when using parallelism, you can use the following - +If the error only occurs when using parallelism, you can use the following example command ```sh mpirun -np 2 xterm -hold -e gdb -ex 'break MPI_Abort' -ex run --args /home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/bin/PHiLiP_2D "-i" "/home/ddong/Codes/PHiLiP_temp/PHiLiP/build_de bug/tests/advection_implicit/2d_advection_implicit_strong.prm" ``` +This launches 2 xterm processes, each of which will launch gdb processes that will run the code and will have a breakpoint when MPI_Abort is encountered. # License From 32eb2299884ff857e8804c13d98471c5611563e6 Mon Sep 17 00:00:00 2001 From: abtin98 Date: Tue, 17 Sep 2019 15:15:51 -0400 Subject: [PATCH 02/10] changed some test parameters --- .../1d_burgers_inviscid_implicit.prm | 4 +++- .../1d_euler_laxfriedrichs_manufactured.prm | 5 +++++ tests/euler_integration/2d_euler_cylinder.prm | 6 ++++++ tests/euler_integration/2d_euler_gaussian_bump.prm | 2 ++ 4 files changed, 16 insertions(+), 1 deletion(-) diff --git a/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm b/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm index 490ea3952..3d2358913 100644 --- a/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm +++ b/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm @@ -5,7 +5,9 @@ set dimension = 1 # The PDE we want to solve. Choices are # . -set pde_type = burgers_inviscid +set pde_type = burgers_inviscid + +set use_weak_form = true subsection ODE solver diff --git a/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm b/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm index 71a33fabb..c92dd89ee 100644 --- a/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm +++ b/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm @@ -7,6 +7,11 @@ set pde_type = euler set conv_num_flux = lax_friedrichs +#set use_split_form = true + +#set use_collocated_nodes = true + +set use_weak_form = true subsection ODE solver diff --git a/tests/euler_integration/2d_euler_cylinder.prm b/tests/euler_integration/2d_euler_cylinder.prm index a97f19d8a..fef66c84a 100644 --- a/tests/euler_integration/2d_euler_cylinder.prm +++ b/tests/euler_integration/2d_euler_cylinder.prm @@ -12,6 +12,12 @@ set pde_type = euler set conv_num_flux = roe +set use_split_form = false + +set use_weak_form = true + +set use_collocated_nodes = false + subsection euler set reference_length = 1.0 set mach_infinity = 0.3 diff --git a/tests/euler_integration/2d_euler_gaussian_bump.prm b/tests/euler_integration/2d_euler_gaussian_bump.prm index f1c2eabb1..74dbad06c 100644 --- a/tests/euler_integration/2d_euler_gaussian_bump.prm +++ b/tests/euler_integration/2d_euler_gaussian_bump.prm @@ -12,6 +12,8 @@ set pde_type = euler set conv_num_flux = roe +set use_split_form = false + subsection euler set reference_length = 1.0 set mach_infinity = 0.3 From f0f3d9e35edd7789d74b5d91a55434a49a2ee129 Mon Sep 17 00:00:00 2001 From: abtin98 Date: Thu, 19 Sep 2019 15:00:50 -0400 Subject: [PATCH 03/10] merged tests with upstream master --- src/parameters/all_parameters.cpp | 10 +- src/testing/CMakeLists.txt | 2 + .../testing}/advection_explicit_periodic.cpp | 184 +++++++------- src/testing/advection_explicit_periodic.h | 26 ++ src/testing/burgers_stability.h | 2 +- ...ler_split_inviscid_taylor_green_vortex.cpp | 238 ++++++++---------- ...euler_split_inviscid_taylor_green_vortex.h | 61 +++++ src/testing/tests.cpp | 6 + .../2D_advection_explicit_periodic.prm | 3 + .../CMakeLists.txt | 59 +---- ...ler_split_inviscid_taylor_green_vortex.prm | 2 + .../CMakeLists.txt | 59 +---- 12 files changed, 315 insertions(+), 337 deletions(-) rename {tests/integration_tests_control_files/advection_explicit_periodic => src/testing}/advection_explicit_periodic.cpp (50%) create mode 100644 src/testing/advection_explicit_periodic.h rename {tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex => src/testing}/euler_split_inviscid_taylor_green_vortex.cpp (55%) create mode 100644 src/testing/euler_split_inviscid_taylor_green_vortex.h diff --git a/src/parameters/all_parameters.cpp b/src/parameters/all_parameters.cpp index 6858c171e..846a77c40 100644 --- a/src/parameters/all_parameters.cpp +++ b/src/parameters/all_parameters.cpp @@ -47,7 +47,9 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm) " euler_vortex | " " euler_entropy_waves | " " numerical_flux_convervation | " - " jacobian_regression "), + " jacobian_regression |" + " advection_periodicity |" + " euler_split_taylor_green"), "The type of test we want to solve. " "Choices are (only run control has been coded up for now)" " ."); + " jacobian_regression |" + " euler_split_taylor_green |" + " advection_periodicity >."); prm.declare_entry("pde_type", "advection", dealii::Patterns::Selection( @@ -111,6 +115,8 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm) else if (test_string == "euler_entropy_waves") { test_type = euler_entropy_waves; } else if (test_string == "numerical_flux_convervation") { test_type = numerical_flux_convervation; } else if (test_string == "jacobian_regression") { test_type = jacobian_regression; } + else if (test_string == "advection_periodicity") {test_type = advection_periodicity; } + else if (test_string == "euler_split_taylor_green") {test_type = euler_split_taylor_green;} const std::string pde_string = prm.get("pde_type"); if (pde_string == "advection") { diff --git a/src/testing/CMakeLists.txt b/src/testing/CMakeLists.txt index e1a22c9c4..ca9a7bd3b 100644 --- a/src/testing/CMakeLists.txt +++ b/src/testing/CMakeLists.txt @@ -6,6 +6,8 @@ set(TESTS_SOURCE euler_cylinder.cpp euler_vortex.cpp euler_entropy_waves.cpp + advection_explicit_periodic.cpp + euler_split_inviscid_taylor_green_vortex.cpp ) foreach(dim RANGE 1 3) diff --git a/tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp b/src/testing/advection_explicit_periodic.cpp similarity index 50% rename from tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp rename to src/testing/advection_explicit_periodic.cpp index 69832e0b4..c3888a047 100644 --- a/tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp +++ b/src/testing/advection_explicit_periodic.cpp @@ -18,43 +18,32 @@ #include "physics/physics.h" #include "dg/dg.h" #include "ode_solver/ode_solver.h" +#include "advection_explicit_periodic.h" #include -template -class AdvectionPeriodic -{ -public: - AdvectionPeriodic() = delete; - AdvectionPeriodic(const PHiLiP::Parameters::AllParameters *const parameters_input); - int run_test(); - -private: - const PHiLiP::Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters - const MPI_Comm mpi_communicator; - dealii::ConditionalOStream pcout; -}; - +namespace PHiLiP { +namespace Tests { template AdvectionPeriodic::AdvectionPeriodic(const PHiLiP::Parameters::AllParameters *const parameters_input) : -all_parameters(parameters_input) +TestsBase::TestsBase(parameters_input) , mpi_communicator(MPI_COMM_WORLD) , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) {} template -int AdvectionPeriodic::run_test() +int AdvectionPeriodic::run_test() const { #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D - dealii::Triangulation grid( - typename dealii::Triangulation::MeshSmoothing( - dealii::Triangulation::smoothing_on_refinement | - dealii::Triangulation::smoothing_on_coarsening)); + dealii::Triangulation grid( + typename dealii::Triangulation::MeshSmoothing( + dealii::Triangulation::smoothing_on_refinement | + dealii::Triangulation::smoothing_on_coarsening)); #else - dealii::parallel::distributed::Triangulation grid( - this->mpi_communicator); + dealii::parallel::distributed::Triangulation grid( + this->mpi_communicator); #endif double left = 0.0; @@ -101,8 +90,8 @@ int AdvectionPeriodic::run_test() constants["pi"] = dealii::numbers::PI; std::string expression = "exp( -( 20*(x-1)*(x-1) + 20*(y-1)*(y-1) ) )";//"sin(pi*x)*sin(pi*y)"; initial_condition.initialize(variables, - expression, - constants); + expression, + constants); dealii::VectorTools::interpolate(dg->dof_handler,initial_condition,dg->solution); // Create ODE solver using the factory and providing the DG object std::shared_ptr> ode_solver = PHiLiP::ODE::ODESolverFactory::create_ODESolver(dg); @@ -115,75 +104,84 @@ int AdvectionPeriodic::run_test() return 0; //need to change } -int main (int argc, char * argv[]) -{ - //parse parameters first - feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan - dealii::deallog.depth_console(99); - dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); - const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); - const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); - dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); - pcout << "Starting program with " << n_mpi << " processors..." << std::endl; - if ((PHILIP_DIM==1) && !(n_mpi==1)) { - std::cout << "********************************************************" << std::endl; - std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl - << "Currently using " << n_mpi << " processors." << std::endl - << "Aborting..." << std::endl; - std::cout << "********************************************************" << std::endl; - std::abort(); - } - int test_error = 1; - try - { - // Declare possible inputs - dealii::ParameterHandler parameter_handler; - PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler); - PHiLiP::Parameters::parse_command_line (argc, argv, parameter_handler); - - // Read inputs from parameter file and set those values in AllParameters object - PHiLiP::Parameters::AllParameters all_parameters; - std::cout << "Reading input..." << std::endl; - all_parameters.parse_parameters (parameter_handler); - - AssertDimension(all_parameters.dimension, PHILIP_DIM); - - std::cout << "Starting program..." << std::endl; - - using namespace PHiLiP; - //const Parameters::AllParameters parameters_input; - AdvectionPeriodic advection_test(&all_parameters); - int i = advection_test.run_test(); - return i; - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl - << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - catch (...) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl - << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - std::cout << "End of program." << std::endl; - return test_error; -} +//int main (int argc, char * argv[]) +//{ +// //parse parameters first +// feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan +// dealii::deallog.depth_console(99); +// dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); +// const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); +// const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); +// dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); +// pcout << "Starting program with " << n_mpi << " processors..." << std::endl; +// if ((PHILIP_DIM==1) && !(n_mpi==1)) { +// std::cout << "********************************************************" << std::endl; +// std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl +// << "Currently using " << n_mpi << " processors." << std::endl +// << "Aborting..." << std::endl; +// std::cout << "********************************************************" << std::endl; +// std::abort(); +// } +// int test_error = 1; +// try +// { +// // Declare possible inputs +// dealii::ParameterHandler parameter_handler; +// PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler); +// PHiLiP::Parameters::parse_command_line (argc, argv, parameter_handler); +// +// // Read inputs from parameter file and set those values in AllParameters object +// PHiLiP::Parameters::AllParameters all_parameters; +// std::cout << "Reading input..." << std::endl; +// all_parameters.parse_parameters (parameter_handler); +// +// AssertDimension(all_parameters.dimension, PHILIP_DIM); +// +// std::cout << "Starting program..." << std::endl; +// +// using namespace PHiLiP; +// //const Parameters::AllParameters parameters_input; +// AdvectionPeriodic advection_test(&all_parameters); +// int i = advection_test.run_test(); +// return i; +// } +// catch (std::exception &exc) +// { +// std::cerr << std::endl << std::endl +// << "----------------------------------------------------" +// << std::endl +// << "Exception on processing: " << std::endl +// << exc.what() << std::endl +// << "Aborting!" << std::endl +// << "----------------------------------------------------" +// << std::endl; +// return 1; +// } +// +// catch (...) +// { +// std::cerr << std::endl +// << std::endl +// << "----------------------------------------------------" +// << std::endl +// << "Unknown exception!" << std::endl +// << "Aborting!" << std::endl +// << "----------------------------------------------------" +// << std::endl; +// return 1; +// } +// std::cout << "End of program." << std::endl; +// return test_error; +//} + +#if PHILIP_DIM==2 + template class AdvectionPeriodic ; +#endif + +} //Tests namespace +} //PHiLiP namespace + + diff --git a/src/testing/advection_explicit_periodic.h b/src/testing/advection_explicit_periodic.h new file mode 100644 index 000000000..c75766878 --- /dev/null +++ b/src/testing/advection_explicit_periodic.h @@ -0,0 +1,26 @@ +#ifndef __ADVECTION_EXPLICIT_PERIODIC_H__ +#define __ADVECTION_EXPLICIT_PERIODIC_H__ + +#include "tests.h" +#include "dg/dg.h" +#include "parameters/all_parameters.h" + +namespace PHiLiP { +namespace Tests { + +template +class AdvectionPeriodic: public TestsBase +{ +public: + AdvectionPeriodic() = delete; + AdvectionPeriodic(const Parameters::AllParameters *const parameters_input); + int run_test () const override; +private: + //const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters + const MPI_Comm mpi_communicator; + dealii::ConditionalOStream pcout; +}; + +} // End of Tests namespace +} // End of PHiLiP namespace +#endif diff --git a/src/testing/burgers_stability.h b/src/testing/burgers_stability.h index f13523f53..4b3bfd703 100644 --- a/src/testing/burgers_stability.h +++ b/src/testing/burgers_stability.h @@ -12,7 +12,7 @@ template class BurgersEnergyStability: public TestsBase { public: - BurgersEnergyStability(const PHiLiP::Parameters::AllParameters *const parameters_input); + BurgersEnergyStability(const Parameters::AllParameters *const parameters_input); int run_test () const override; private: double compute_energy(std::shared_ptr < PHiLiP::DGBase > &dg) const; diff --git a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp b/src/testing/euler_split_inviscid_taylor_green_vortex.cpp similarity index 55% rename from tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp rename to src/testing/euler_split_inviscid_taylor_green_vortex.cpp index 9d104faee..5461da01e 100644 --- a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp +++ b/src/testing/euler_split_inviscid_taylor_green_vortex.cpp @@ -1,77 +1,29 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - - -#include "parameters/all_parameters.h" -#include "parameters/parameters.h" -#include "numerical_flux/numerical_flux.h" -#include "physics/physics_factory.h" -#include "physics/physics.h" -#include "dg/dg.h" -#include "ode_solver/ode_solver.h" - -#include - -//using PDEType = PHiLiP::Parameters::AllParameters::PartialDifferentialEquation; -//using ConvType = PHiLiP::Parameters::AllParameters::ConvectiveNumericalFlux; -//using DissType = PHiLiP::Parameters::AllParameters::DissipativeNumericalFlux; -// -// -//const double TOLERANCE = 1E-12; - - -template -class EulerTaylorGreen -{ -public: - EulerTaylorGreen() = delete; - EulerTaylorGreen(const PHiLiP::Parameters::AllParameters *const parameters_input); - int run_test(); - - - +#include "euler_split_inviscid_taylor_green_vortex.h" -private: - double compute_kinetic_energy(std::shared_ptr < PHiLiP::DGBase > &dg, unsigned int poly_degree); - double compute_quadrature_kinetic_energy(std::array soln_at_q); - const PHiLiP::Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters - const MPI_Comm mpi_communicator; - dealii::ConditionalOStream pcout; -}; +namespace PHiLiP { +namespace Tests { template -EulerTaylorGreen::EulerTaylorGreen(const PHiLiP::Parameters::AllParameters *const parameters_input) +EulerTaylorGreen::EulerTaylorGreen(const Parameters::AllParameters *const parameters_input) : -all_parameters(parameters_input) +TestsBase::TestsBase(parameters_input) , mpi_communicator(MPI_COMM_WORLD) , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) {} -template -double EulerTaylorGreen::compute_quadrature_kinetic_energy(std::array soln_at_q) -{ - const double density = soln_at_q[0]; - - const double quad_kin_energ = 0.5*(soln_at_q[1]*soln_at_q[1] + - soln_at_q[2]*soln_at_q[2] + - soln_at_q[3]*soln_at_q[3] )/density; - return quad_kin_energ; -} +//template +//double EulerTaylorGreen::compute_quadrature_kinetic_energy(std::array soln_at_q) const +//{ +// const double density = soln_at_q[0]; +// +// const double quad_kin_energ = 0.5*(soln_at_q[1]*soln_at_q[1] + +// soln_at_q[2]*soln_at_q[2] + +// soln_at_q[3]*soln_at_q[3] )/density; +// return quad_kin_energ; +//} template -double EulerTaylorGreen::compute_kinetic_energy(std::shared_ptr < PHiLiP::DGBase > &dg, unsigned int poly_degree) +double EulerTaylorGreen::compute_kinetic_energy(std::shared_ptr < DGBase > &dg, unsigned int poly_degree) const { // Overintegrate the error to make sure there is not integration error in the error estimate int overintegrate = 10 ;//10; @@ -114,7 +66,14 @@ double EulerTaylorGreen::compute_kinetic_energy(std::shared_ptr < P const unsigned int istate = fe_values_extra.get_fe().system_to_component_index(idof).first; soln_at_q[istate] += dg->solution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate); } - const double quadrature_kinetic_energy = compute_quadrature_kinetic_energy(soln_at_q); + + const double density = soln_at_q[0]; + + const double quadrature_kinetic_energy = 0.5*(soln_at_q[1]*soln_at_q[1] + + soln_at_q[2]*soln_at_q[2] + + soln_at_q[3]*soln_at_q[3] )/density; + + //const double quadrature_kinetic_energy = compute_quadrature_kinetic_energy(soln_at_q); total_kinetic_energy += quadrature_kinetic_energy * fe_values_extra.JxW(iquad); } @@ -123,7 +82,7 @@ double EulerTaylorGreen::compute_kinetic_energy(std::shared_ptr < P } template -int EulerTaylorGreen::run_test() +int EulerTaylorGreen::run_test() const { //dealii::Triangulation grid; dealii::parallel::distributed::Triangulation grid(this->mpi_communicator); @@ -143,7 +102,7 @@ int EulerTaylorGreen::run_test() grid.refine_global(n_refinements); - std::shared_ptr < PHiLiP::DGBase > dg = PHiLiP::DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree); dg->set_triangulation(&grid); dg->allocate_system (); @@ -172,7 +131,7 @@ int EulerTaylorGreen::run_test() // Create ODE solver using the factory and providing the DG object std::cout << "creating ODE solver" << std::endl; - std::shared_ptr> ode_solver = PHiLiP::ODE::ODESolverFactory::create_ODESolver(dg); + std::shared_ptr> ode_solver = ODE::ODESolverFactory::create_ODESolver(dg); std::cout << "ODE solver successfully created" << std::endl; double finalTime = 14.; double dt = all_parameters->ode_solver_param.initial_time_step; @@ -213,77 +172,82 @@ int EulerTaylorGreen::run_test() return 0; //need to change } -int main (int argc, char * argv[]) -{ - //parse parameters first - feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan - dealii::deallog.depth_console(99); - dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); - const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); - const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); - dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); - pcout << "Starting program with " << n_mpi << " processors..." << std::endl; - if ((PHILIP_DIM==1) && !(n_mpi==1)) { - std::cout << "********************************************************" << std::endl; - std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl - << "Currently using " << n_mpi << " processors." << std::endl - << "Aborting..." << std::endl; - std::cout << "********************************************************" << std::endl; - std::abort(); - } - int test_error = 1; - try - { - // Declare possible inputs - dealii::ParameterHandler parameter_handler; - PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler); - PHiLiP::Parameters::parse_command_line (argc, argv, parameter_handler); - - // Read inputs from parameter file and set those values in AllParameters object - PHiLiP::Parameters::AllParameters all_parameters; - std::cout << "Reading input..." << std::endl; - all_parameters.parse_parameters (parameter_handler); - - AssertDimension(all_parameters.dimension, PHILIP_DIM); - - std::cout << "Starting program..." << std::endl; - - using namespace PHiLiP; - //const Parameters::AllParameters parameters_input; - EulerTaylorGreen euler_test(&all_parameters); - int i = euler_test.run_test(); - return i; - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl - << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } +//int main (int argc, char * argv[]) +//{ +// //parse parameters first +// feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan +// dealii::deallog.depth_console(99); +// dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); +// const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); +// const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); +// dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); +// pcout << "Starting program with " << n_mpi << " processors..." << std::endl; +// if ((PHILIP_DIM==1) && !(n_mpi==1)) { +// std::cout << "********************************************************" << std::endl; +// std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl +// << "Currently using " << n_mpi << " processors." << std::endl +// << "Aborting..." << std::endl; +// std::cout << "********************************************************" << std::endl; +// std::abort(); +// } +// int test_error = 1; +// try +// { +// // Declare possible inputs +// dealii::ParameterHandler parameter_handler; +// Parameters::AllParameters::declare_parameters (parameter_handler); +// Parameters::parse_command_line (argc, argv, parameter_handler); +// +// // Read inputs from parameter file and set those values in AllParameters object +// Parameters::AllParameters all_parameters; +// std::cout << "Reading input..." << std::endl; +// all_parameters.parse_parameters (parameter_handler); +// +// AssertDimension(all_parameters.dimension, PHILIP_DIM); +// +// std::cout << "Starting program..." << std::endl; +// +// using namespace PHiLiP; +// //const Parameters::AllParameters parameters_input; +// EulerTaylorGreen euler_test(&all_parameters); +// int i = euler_test.run_test(); +// return i; +// } +// catch (std::exception &exc) +// { +// std::cerr << std::endl << std::endl +// << "----------------------------------------------------" +// << std::endl +// << "Exception on processing: " << std::endl +// << exc.what() << std::endl +// << "Aborting!" << std::endl +// << "----------------------------------------------------" +// << std::endl; +// return 1; +// } +// +// catch (...) +// { +// std::cerr << std::endl +// << std::endl +// << "----------------------------------------------------" +// << std::endl +// << "Unknown exception!" << std::endl +// << "Aborting!" << std::endl +// << "----------------------------------------------------" +// << std::endl; +// return 1; +// } +// std::cout << "End of program." << std::endl; +// return test_error; +//} + +#if PHILIP_DIM==3 + template class EulerTaylorGreen ; +#endif - catch (...) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl - << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - std::cout << "End of program." << std::endl; - return test_error; } - - +} diff --git a/src/testing/euler_split_inviscid_taylor_green_vortex.h b/src/testing/euler_split_inviscid_taylor_green_vortex.h new file mode 100644 index 000000000..58a3d9314 --- /dev/null +++ b/src/testing/euler_split_inviscid_taylor_green_vortex.h @@ -0,0 +1,61 @@ +#ifndef __EULER_SPLIT_TAYLOR_GREEN_H__ +#define __EULER_SPLIT_TAYLOR_GREEN_H__ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include "tests.h" + + +#include "parameters/all_parameters.h" +#include "parameters/parameters.h" +#include "numerical_flux/numerical_flux.h" +#include "physics/physics_factory.h" +#include "physics/physics.h" +#include "dg/dg.h" +#include "ode_solver/ode_solver.h" + +#include + +//using PDEType = PHiLiP::Parameters::AllParameters::PartialDifferentialEquation; +//using ConvType = PHiLiP::Parameters::AllParameters::ConvectiveNumericalFlux; +//using DissType = PHiLiP::Parameters::AllParameters::DissipativeNumericalFlux; +// +// +//const double TOLERANCE = 1E-12; +namespace PHiLiP { +namespace Tests { + +template +class EulerTaylorGreen : public TestsBase +{ +public: + EulerTaylorGreen() = delete; + EulerTaylorGreen(const Parameters::AllParameters *const parameters_input); + int run_test() const override; + +private: + double compute_kinetic_energy(std::shared_ptr < DGBase > &dg, unsigned int poly_degree) const; + //double compute_quadrature_kinetic_energy(std::array soln_at_q) const ; + //const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters + const MPI_Comm mpi_communicator; + dealii::ConditionalOStream pcout; +}; + + +} //Tests +} //PHiLiP + +#endif diff --git a/src/testing/tests.cpp b/src/testing/tests.cpp index c92f6f03a..2cd8e5142 100644 --- a/src/testing/tests.cpp +++ b/src/testing/tests.cpp @@ -11,6 +11,8 @@ #include "euler_cylinder.h" #include "euler_vortex.h" #include "euler_entropy_waves.h" +#include "advection_explicit_periodic.h" +#include "euler_split_inviscid_taylor_green_vortex.h" namespace PHiLiP { namespace Tests { @@ -57,6 +59,8 @@ ::select_test(const AllParam *const parameters_input) { return std::make_unique>(parameters_input); } else if(test_type == Test_enum::burgers_energy_stability) { if constexpr (dim==1 && nstate==1) return std::make_unique>(parameters_input); + } else if (test_type == Test_enum::advection_periodicity){ + if constexpr (dim == 2 && nstate == 1) return std::make_unique> (parameters_input); } else if(test_type == Test_enum::euler_gaussian_bump) { if constexpr (dim==2 && nstate==dim+2) return std::make_unique>(parameters_input); } else if(test_type == Test_enum::euler_cylinder) { @@ -65,6 +69,8 @@ ::select_test(const AllParam *const parameters_input) { if constexpr (dim==2 && nstate==dim+2) return std::make_unique>(parameters_input); } else if(test_type == Test_enum::euler_entropy_waves) { if constexpr (dim>=2 && nstate==PHILIP_DIM+2) return std::make_unique>(parameters_input); + } else if(test_type == Test_enum::euler_split_taylor_green) { + if constexpr (dim==3 && nstate == dim+2) return std::make_unique>(parameters_input); } else { std::cout << "Invalid test." << std::endl; } diff --git a/tests/integration_tests_control_files/advection_explicit_periodic/2D_advection_explicit_periodic.prm b/tests/integration_tests_control_files/advection_explicit_periodic/2D_advection_explicit_periodic.prm index 8b4798e6d..7b3a7a4b7 100644 --- a/tests/integration_tests_control_files/advection_explicit_periodic/2D_advection_explicit_periodic.prm +++ b/tests/integration_tests_control_files/advection_explicit_periodic/2D_advection_explicit_periodic.prm @@ -1,4 +1,7 @@ # ------------------- + +set test_type = advection_periodicity + # Number of dimensions set dimension = 2 diff --git a/tests/integration_tests_control_files/advection_explicit_periodic/CMakeLists.txt b/tests/integration_tests_control_files/advection_explicit_periodic/CMakeLists.txt index 78ee6c466..35c2766ec 100644 --- a/tests/integration_tests_control_files/advection_explicit_periodic/CMakeLists.txt +++ b/tests/integration_tests_control_files/advection_explicit_periodic/CMakeLists.txt @@ -1,52 +1,7 @@ -set(TEST_SRC - advection_explicit_periodic.cpp - ) - -foreach(dim RANGE 2 2) - - - # Output executable - string(CONCAT TEST_TARGET ${dim}D_advection_explicit_periodic) - message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n") - add_executable(${TEST_TARGET} ${TEST_SRC}) - # Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code - target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim}) - - # Compile this executable when 'make unit_tests' - add_dependencies(unit_tests ${TEST_TARGET}) - add_dependencies(${dim}D ${TEST_TARGET}) - - # Library dependency - set(ParameterLib ParametersLibrary) - string(CONCAT PhysicsLib Physics_${dim}D) - string(CONCAT NumericalFluxLib NumericalFlux_${dim}D) - string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) - string(CONCAT ODESolverLib ODESolver_${dim}D) - - target_link_libraries(${TEST_TARGET} ${ParameterLib}) - target_link_libraries(${TEST_TARGET} ${PhysicsLib}) - target_link_libraries(${TEST_TARGET} ${NumericalFluxLib}) - target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) - target_link_libraries(${TEST_TARGET} ${ODESolverLib}) - - #target_link_libraries(${TEST_TARGET} m mpi mpi_cxx) - - # Setup target with deal.II - DEAL_II_SETUP_TARGET(${TEST_TARGET}) - - configure_file(2D_advection_explicit_periodic.prm 2D_advection_explicit_periodic.prm COPYONLY) - add_test( - NAME ${TEST_TARGET} - COMMAND mpirun -np ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} -i ${CMAKE_CURRENT_BINARY_DIR}/2D_advection_explicit_periodic.prm - WORKING_DIRECTORY ${TEST_OUTPUT_DIR} - ) - - unset(dim) - unset(TEST_TARGET) - unset(PhysicsLib) - unset(NumericalFluxLib) - unset(ParameterLib) - unset(DiscontinuousGalerkinLib) - unset(ODESolverLib) - -endforeach() +set(TEST_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}) +configure_file(2D_advection_explicit_periodic.prm 2D_advection_explicit_periodic.prm COPYONLY) +add_test( + NAME MPI_2D_ADVECTION_EXPLICIT_PERIODIC +COMMAND mpirun -n ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/2D_advection_explicit_periodic.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) diff --git a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm b/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm index a3946e579..bac24f4c9 100644 --- a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm +++ b/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm @@ -2,6 +2,8 @@ # Number of dimensions set dimension = 3 +set test_type = euler_split_taylor_green + set use_weak_form = false set use_collocated_nodes = true diff --git a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt b/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt index 215be90b0..f4c9930c5 100644 --- a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt +++ b/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt @@ -1,52 +1,7 @@ -set(TEST_SRC - euler_split_inviscid_taylor_green_vortex.cpp - ) - -foreach(dim RANGE 3 3) - - - # Output executable - string(CONCAT TEST_TARGET ${dim}D_euler_split_inviscid_taylor_green_vortex) - message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n") - add_executable(${TEST_TARGET} ${TEST_SRC}) - # Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code - target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim}) - - # Compile this executable when 'make unit_tests' - add_dependencies(unit_tests ${TEST_TARGET}) - add_dependencies(${dim}D ${TEST_TARGET}) - - # Library dependency - set(ParameterLib ParametersLibrary) - string(CONCAT PhysicsLib Physics_${dim}D) - string(CONCAT NumericalFluxLib NumericalFlux_${dim}D) - string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) - string(CONCAT ODESolverLib ODESolver_${dim}D) - - target_link_libraries(${TEST_TARGET} ${ParameterLib}) - target_link_libraries(${TEST_TARGET} ${PhysicsLib}) - target_link_libraries(${TEST_TARGET} ${NumericalFluxLib}) - target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) - target_link_libraries(${TEST_TARGET} ${ODESolverLib}) - - #target_link_libraries(${TEST_TARGET} m mpi mpi_cxx) - - # Setup target with deal.II - DEAL_II_SETUP_TARGET(${TEST_TARGET}) - - configure_file(3D_euler_split_inviscid_taylor_green_vortex.prm 3D_euler_split_inviscid_taylor_green_vortex.prm COPYONLY) - add_test( - NAME ${TEST_TARGET} - COMMAND mpirun -np 1 ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} -i ${CMAKE_CURRENT_BINARY_DIR}/3D_euler_split_inviscid_taylor_green_vortex.prm - WORKING_DIRECTORY ${TEST_OUTPUT_DIR} - ) - - unset(dim) - unset(TEST_TARGET) - unset(PhysicsLib) - unset(NumericalFluxLib) - unset(ParameterLib) - unset(DiscontinuousGalerkinLib) - unset(ODESolverLib) - -endforeach() +set(TEST_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}) +configure_file(3D_euler_split_inviscid_taylor_green_vortex.prm 3D_euler_split_inviscid_taylor_green_vortex.prm COPYONLY) +add_test( + NAME MPI_3D_EULER_SPLIT_TAYLOR_GREEN + COMMAND mpirun -n ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/3D_euler_split_inviscid_taylor_green_vortex.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) From 2681a89d30bb6c3bb8bc705dae18c3cdb6fed99c Mon Sep 17 00:00:00 2001 From: abtin98 Date: Mon, 23 Sep 2019 17:55:29 -0400 Subject: [PATCH 04/10] fixed collocation bug --- src/dg/dg.cpp | 60 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp index f4ad0077d..b260e73f3 100644 --- a/src/dg/dg.cpp +++ b/src/dg/dg.cpp @@ -129,6 +129,66 @@ DGBase::create_collection_tuple(const unsigned int max_degree, const i dealii::hp::QCollection<1> oned_quad_coll; dealii::hp::FECollection fe_coll_lagr; + + // for p=0, we use a p=1 FE for collocation, since there's no p=0 quadrature for Gauss Lobatto + if (parameters_input->use_collocated_nodes==true) + { + int degree = 1; + //const dealii::MappingQ mapping(degree, true); + //const dealii::MappingQ mapping(degree+1, true); + const dealii::MappingManifold mapping; + mapping_coll.push_back(mapping); + + const dealii::FE_DGQ fe_dg(degree); + const dealii::FESystem fe_system(fe_dg, nstate); + fe_coll.push_back (fe_system); + + // + + dealii::Quadrature<1> oned_quad(degree+1); + dealii::Quadrature volume_quad(degree+1); + dealii::Quadrature face_quad(degree+1); //removed const + + if (parameters_input->use_collocated_nodes) + { + dealii::QGaussLobatto<1> oned_quad_Gauss_Lobatto (degree+1); + dealii::QGaussLobatto vol_quad_Gauss_Lobatto (degree+1); + oned_quad = oned_quad_Gauss_Lobatto; + volume_quad = vol_quad_Gauss_Lobatto; + + if(dim == 1) + { + dealii::QGauss face_quad_Gauss_Legendre (degree+1); + face_quad = face_quad_Gauss_Legendre; + } + else + { + dealii::QGaussLobatto face_quad_Gauss_Lobatto (degree+1); + face_quad = face_quad_Gauss_Lobatto; + } + + + } + else + { + dealii::QGauss<1> oned_quad_Gauss_Legendre (degree+1); + dealii::QGauss vol_quad_Gauss_Legendre (degree+1); + dealii::QGauss face_quad_Gauss_Legendre (degree+1); + oned_quad = oned_quad_Gauss_Legendre; + volume_quad = vol_quad_Gauss_Legendre; + face_quad = face_quad_Gauss_Legendre; + } + // + + + volume_quad_coll.push_back (volume_quad); + face_quad_coll.push_back (face_quad); + oned_quad_coll.push_back (oned_quad); + + dealii::FE_DGQArbitraryNodes lagrange_poly(oned_quad); + fe_coll_lagr.push_back (lagrange_poly); + } + int minimum_degree = (parameters_input->use_collocated_nodes==true) ? 1 : 0; for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) { //const dealii::MappingQ mapping(degree, true); From e6ebeaaf0f32d52338ff5a75ebafe01356d5a59e Mon Sep 17 00:00:00 2001 From: abtin98 Date: Mon, 23 Sep 2019 17:55:54 -0400 Subject: [PATCH 05/10] added some parts of MHD --- src/physics/CMakeLists.txt | 2 +- src/physics/mhd.cpp | 678 ++++++++++++++++---------------- src/physics/mhd.h | 40 +- src/physics/physics.cpp | 2 + src/physics/physics_factory.cpp | 10 +- 5 files changed, 369 insertions(+), 363 deletions(-) diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt index 91f987ff4..b3a588035 100644 --- a/src/physics/CMakeLists.txt +++ b/src/physics/CMakeLists.txt @@ -5,7 +5,7 @@ SET(SOURCE burgers.cpp euler.cpp manufactured_solution.cpp - # mhd.cpp + mhd.cpp ) foreach(dim RANGE 1 3) diff --git a/src/physics/mhd.cpp b/src/physics/mhd.cpp index 6d4b404a2..c1decd7e8 100644 --- a/src/physics/mhd.cpp +++ b/src/physics/mhd.cpp @@ -27,6 +27,7 @@ ::source_term ( return source_term; } +//incomplete template inline std::array MHD ::convert_conservative_to_primitive ( const std::array &conservative_soln ) const @@ -45,6 +46,7 @@ ::convert_conservative_to_primitive ( const std::array &conservativ return primitive_soln; } +//incomplete template inline std::array MHD ::convert_primitive_to_conservative ( const std::array &primitive_soln ) const @@ -63,13 +65,6 @@ ::convert_primitive_to_conservative ( const std::array &primitive_s return conservative_soln; } -//template -//inline dealii::Tensor<1,dim,double> MHD::compute_velocities_inf() const -//{ -// dealii::Tensor<1,dim,double> velocities; -// return velocities; -//} - template inline dealii::Tensor<1,dim,real> MHD ::compute_velocities ( const std::array &conservative_soln ) const @@ -98,6 +93,8 @@ ::extract_velocities_from_primitive ( const std::array &primitive_s return velocities; } + +//incomplete template inline real MHD ::compute_total_energy ( const std::array &primitive_soln ) const @@ -111,6 +108,7 @@ ::compute_total_energy ( const std::array &primitive_soln ) const return tot_energy; } +//incomplete template inline real MHD ::compute_entropy_measure ( const std::array &conservative_soln ) const @@ -121,7 +119,7 @@ ::compute_entropy_measure ( const std::array &conservative_soln ) c return entropy_measure; } - +//incomplete template inline real MHD ::compute_specific_enthalpy ( const std::array &conservative_soln, const real pressure ) const @@ -132,6 +130,7 @@ ::compute_specific_enthalpy ( const std::array &conservative_soln, return specific_enthalpy; } + template inline real MHD ::compute_dimensional_temperature ( const std::array &primitive_soln ) const @@ -147,7 +146,7 @@ inline real MHD ::compute_temperature ( const std::array &primitive_soln ) const { const real dimensional_temperature = compute_dimensional_temperature(primitive_soln); - const real temperature = dimensional_temperature * mach_inf_sqr; + const real temperature = dimensional_temperature /** mach_inf_sqr*/; return temperature; } @@ -155,14 +154,14 @@ template inline real MHD ::compute_density_from_pressure_temperature ( const real pressure, const real temperature ) const { - const real density = gam*pressure/temperature * mach_inf_sqr; + const real density = gam*pressure/temperature /** mach_inf_sqr*/; return density; } template inline real MHD ::compute_temperature_from_density_pressure ( const real density, const real pressure ) const { - const real temperature = gam*pressure/density * mach_inf_sqr; + const real temperature = gam*pressure/density /* mach_inf_sqr*/; return temperature; } @@ -239,9 +238,10 @@ template inline real MHD ::compute_magnetic_energy (const std::array &conservative_soln) const { - const real magnetic_energy = 0; + real magnetic_energy = 0; for (int i = 1; i <= 3; ++i) magnetic_energy += 1./2. * (conservative_soln[nstate - i] * conservative_soln[nstate - i] ); + return magnetic_energy; } @@ -463,335 +463,335 @@ ::dissipative_flux ( return diss_flux; } -template -void MHD -::boundary_face_values ( - const int boundary_type, - const dealii::Point &pos, - const dealii::Tensor<1,dim,real> &normal_int, - const std::array &soln_int, - const std::array,nstate> &soln_grad_int, - std::array &soln_bc, - std::array,nstate> &soln_grad_bc) const -{ - // NEED TO PROVIDE AS INPUT ************************************** - const real total_inlet_pressure = pressure_inf*pow(1.0+0.5*gamm1*mach_inf_sqr, gam/gamm1); - const real total_inlet_temperature = temperature_inf*pow(total_inlet_pressure/pressure_inf, gamm1/gam); - - if (boundary_type == 1000) { - // Manufactured solution - std::array conservative_boundary_values; - std::array,nstate> boundary_gradients; - for (int s=0; smanufactured_solution_function.value (pos, s); - boundary_gradients[s] = this->manufactured_solution_function.gradient (pos, s); - } - std::array primitive_boundary_values = convert_conservative_to_primitive(conservative_boundary_values); - for (int istate=0; istate characteristic_dot_n = convective_eigenvalues(conservative_boundary_values, normal_int); - const bool inflow = (characteristic_dot_n[istate] <= 0.); - - if (inflow) { // Dirichlet boundary condition - - soln_bc[istate] = conservative_boundary_values[istate]; - soln_grad_bc[istate] = soln_grad_int[istate]; - - // Only set the pressure and velocity - // primitive_boundary_values[0] = soln_int[0];; - // for(int d=0;d primitive_interior_values = convert_conservative_to_primitive(soln_int); - - // Copy density and pressure - std::array primitive_boundary_values; - primitive_boundary_values[0] = primitive_interior_values[0]; - primitive_boundary_values[nstate-1] = primitive_interior_values[nstate-1]; - - const dealii::Tensor<1,dim,real> surface_normal = -normal_int; - const dealii::Tensor<1,dim,real> velocities_int = extract_velocities_from_primitive(primitive_interior_values); - const dealii::Tensor<1,dim,real> velocities_bc = velocities_int - 2.0*(velocities_int*surface_normal)*surface_normal; - for (int d=0; d primitive_interior_values = convert_conservative_to_primitive(soln_int); - const real pressure_int = primitive_interior_values[nstate-1]; - - const real radicant = 1.0+0.5*gamm1*mach_inf_sqr; - const real pressure_inlet = total_inlet_pressure * pow(radicant, -gam/gamm1); - const real pressure_bc = (mach_int >= 1) ? pressure_int : back_pressure*pressure_inlet; - const real temperature_int = compute_temperature(primitive_interior_values); - - // Assign primitive boundary values - std::array primitive_boundary_values; - primitive_boundary_values[0] = compute_density_from_pressure_temperature(pressure_bc, temperature_int); - for (int d=0;d 1.0) { - soln_bc = soln_int; - } - - } else if (boundary_type == 1003) { - // Inflow - // Carlson 2011, sec. 2.2 & sec 2.9 - - const std::array primitive_interior_values = convert_conservative_to_primitive(soln_int); - - const dealii::Tensor<1,dim,real> normal = -normal_int; - - const real density_i = primitive_interior_values[0]; - const dealii::Tensor<1,dim,real> velocities_i = extract_velocities_from_primitive(primitive_interior_values); - const real pressure_i = primitive_interior_values[nstate-1]; - - const real normal_vel_i = velocities_i*normal; - const real sound_i = compute_sound(soln_int); - //const real mach_i = std::abs(normal_vel_i)/sound_i; - - //const dealii::Tensor<1,dim,real> velocities_o = velocities_inf; - //const real normal_vel_o = velocities_o*normal; - //const real sound_o = sound_inf; - //const real mach_o = mach_inf; - - if(mach_inf < 1.0) { - //std::cout << "Subsonic inflow, mach=" << mach_i << std::endl; - // Subsonic inflow, sec 2.7 - - // Want to solve for c_b (sound_bc), to then solve for U (velocity_magnitude_bc) and M_b (mach_bc) - // Eq. 37 - const real riemann_pos = normal_vel_i + 2.0*sound_i/gamm1; - // Could evaluate enthalpy from primitive like eq.36, but easier to use the following - const real specific_total_energy = soln_int[nstate-1]/density_i; - const real specific_total_enthalpy = specific_total_energy + pressure_i/density_i; - // Eq. 43 - const real a = 1.0+2.0/gamm1; - const real b = -2.0*riemann_pos; - const real c = 0.5*gamm1 * (riemann_pos*riemann_pos - 2.0*specific_total_enthalpy); - // Eq. 42 - const real term1 = -0.5*b/a; - const real term2= 0.5*sqrt(b*b-4.0*a*c)/a; - const real sound_bc1 = term1 + term2; - const real sound_bc2 = term1 - term2; - // Eq. 44 - const real sound_bc = std::max(sound_bc1, sound_bc2); - // Eq. 45 - //const real velocity_magnitude_bc = 2.0*sound_bc/gamm1 - riemann_pos; - const real velocity_magnitude_bc = riemann_pos - 2.0*sound_bc/gamm1; - const real mach_bc = velocity_magnitude_bc/sound_bc; - // Eq. 46 - const real radicant = 1.0+0.5*gamm1*mach_bc*mach_bc; - const real pressure_bc = total_inlet_pressure * pow(radicant, -gam/gamm1); - const real temperature_bc = total_inlet_temperature * pow(radicant, -1.0); - //std::cout << " pressure_bc " << pressure_bc << "pressure_inf" << pressure_inf << std::endl; - //std::cout << " temperature_bc " << temperature_bc << "temperature_inf" << temperature_inf << std::endl; - // - - const real density_bc = compute_density_from_pressure_temperature(pressure_bc, temperature_bc); - std::array primitive_boundary_values; - primitive_boundary_values[0] = density_bc; - for (int d=0;d primitive_boundary_values; - primitive_boundary_values[0] = density_bc; - for (int d=0;d -dealii::Vector MHD::post_compute_derived_quantities_vector ( - const dealii::Vector &uh, - const std::vector > &duh, - const std::vector > &dduh, - const dealii::Tensor<1,dim> &normals, - const dealii::Point &evaluation_points) const -{ - std::vector names = post_get_names (); - dealii::Vector computed_quantities = PhysicsBase::post_compute_derived_quantities_vector ( uh, duh, dduh, normals, evaluation_points); - unsigned int current_data_index = computed_quantities.size() - 1; - computed_quantities.grow_or_shrink(names.size()); - if constexpr (std::is_same::value) { - - std::array conservative_soln; - for (unsigned int s=0; s primitive_soln = convert_conservative_to_primitive(conservative_soln); - - // Density - computed_quantities(++current_data_index) = primitive_soln[0]; - // Velocities - for (unsigned int d=0; d -std::vector MHD -::post_get_data_component_interpretation () const -{ - namespace DCI = dealii::DataComponentInterpretation; - std::vector interpretation = PhysicsBase::post_get_data_component_interpretation (); // state variables - interpretation.push_back (DCI::component_is_scalar); // Density - for (unsigned int d=0; d names = post_get_names(); - if (names.size() != interpretation.size()) { - std::cout << "Number of DataComponentInterpretation is not the same as number of names for output file" << std::endl; - } - return interpretation; -} - - -template -std::vector MHD ::post_get_names () const -{ - std::vector names = PhysicsBase::post_get_names (); - names.push_back ("density"); - for (unsigned int d=0; d -dealii::UpdateFlags MHD -::post_get_needed_update_flags () const -{ - //return update_values | update_gradients; - return dealii::update_values; -} +//template +//void MHD +//::boundary_face_values ( +// const int boundary_type, +// const dealii::Point &pos, +// const dealii::Tensor<1,dim,real> &normal_int, +// const std::array &soln_int, +// const std::array,nstate> &soln_grad_int, +// std::array &soln_bc, +// std::array,nstate> &soln_grad_bc) const +//{ +// // NEED TO PROVIDE AS INPUT ************************************** +// const real total_inlet_pressure = pressure_inf*pow(1.0+0.5*gamm1*mach_inf_sqr, gam/gamm1); +// const real total_inlet_temperature = temperature_inf*pow(total_inlet_pressure/pressure_inf, gamm1/gam); +// +// if (boundary_type == 1000) { +// // Manufactured solution +// std::array conservative_boundary_values; +// std::array,nstate> boundary_gradients; +// for (int s=0; smanufactured_solution_function.value (pos, s); +// boundary_gradients[s] = this->manufactured_solution_function.gradient (pos, s); +// } +// std::array primitive_boundary_values = convert_conservative_to_primitive(conservative_boundary_values); +// for (int istate=0; istate characteristic_dot_n = convective_eigenvalues(conservative_boundary_values, normal_int); +// const bool inflow = (characteristic_dot_n[istate] <= 0.); +// +// if (inflow) { // Dirichlet boundary condition +// +// soln_bc[istate] = conservative_boundary_values[istate]; +// soln_grad_bc[istate] = soln_grad_int[istate]; +// +// // Only set the pressure and velocity +// // primitive_boundary_values[0] = soln_int[0];; +// // for(int d=0;d primitive_interior_values = convert_conservative_to_primitive(soln_int); +// +// // Copy density and pressure +// std::array primitive_boundary_values; +// primitive_boundary_values[0] = primitive_interior_values[0]; +// primitive_boundary_values[nstate-1] = primitive_interior_values[nstate-1]; +// +// const dealii::Tensor<1,dim,real> surface_normal = -normal_int; +// const dealii::Tensor<1,dim,real> velocities_int = extract_velocities_from_primitive(primitive_interior_values); +// const dealii::Tensor<1,dim,real> velocities_bc = velocities_int - 2.0*(velocities_int*surface_normal)*surface_normal; +// for (int d=0; d primitive_interior_values = convert_conservative_to_primitive(soln_int); +// const real pressure_int = primitive_interior_values[nstate-1]; +// +// const real radicant = 1.0+0.5*gamm1*mach_inf_sqr; +// const real pressure_inlet = total_inlet_pressure * pow(radicant, -gam/gamm1); +// const real pressure_bc = (mach_int >= 1) ? pressure_int : back_pressure*pressure_inlet; +// const real temperature_int = compute_temperature(primitive_interior_values); +// +// // Assign primitive boundary values +// std::array primitive_boundary_values; +// primitive_boundary_values[0] = compute_density_from_pressure_temperature(pressure_bc, temperature_int); +// for (int d=0;d 1.0) { +// soln_bc = soln_int; +// } +// +// } else if (boundary_type == 1003) { +// // Inflow +// // Carlson 2011, sec. 2.2 & sec 2.9 +// +// const std::array primitive_interior_values = convert_conservative_to_primitive(soln_int); +// +// const dealii::Tensor<1,dim,real> normal = -normal_int; +// +// const real density_i = primitive_interior_values[0]; +// const dealii::Tensor<1,dim,real> velocities_i = extract_velocities_from_primitive(primitive_interior_values); +// const real pressure_i = primitive_interior_values[nstate-1]; +// +// const real normal_vel_i = velocities_i*normal; +// const real sound_i = compute_sound(soln_int); +// //const real mach_i = std::abs(normal_vel_i)/sound_i; +// +// //const dealii::Tensor<1,dim,real> velocities_o = velocities_inf; +// //const real normal_vel_o = velocities_o*normal; +// //const real sound_o = sound_inf; +// //const real mach_o = mach_inf; +// +// if(mach_inf < 1.0) { +// //std::cout << "Subsonic inflow, mach=" << mach_i << std::endl; +// // Subsonic inflow, sec 2.7 +// +// // Want to solve for c_b (sound_bc), to then solve for U (velocity_magnitude_bc) and M_b (mach_bc) +// // Eq. 37 +// const real riemann_pos = normal_vel_i + 2.0*sound_i/gamm1; +// // Could evaluate enthalpy from primitive like eq.36, but easier to use the following +// const real specific_total_energy = soln_int[nstate-1]/density_i; +// const real specific_total_enthalpy = specific_total_energy + pressure_i/density_i; +// // Eq. 43 +// const real a = 1.0+2.0/gamm1; +// const real b = -2.0*riemann_pos; +// const real c = 0.5*gamm1 * (riemann_pos*riemann_pos - 2.0*specific_total_enthalpy); +// // Eq. 42 +// const real term1 = -0.5*b/a; +// const real term2= 0.5*sqrt(b*b-4.0*a*c)/a; +// const real sound_bc1 = term1 + term2; +// const real sound_bc2 = term1 - term2; +// // Eq. 44 +// const real sound_bc = std::max(sound_bc1, sound_bc2); +// // Eq. 45 +// //const real velocity_magnitude_bc = 2.0*sound_bc/gamm1 - riemann_pos; +// const real velocity_magnitude_bc = riemann_pos - 2.0*sound_bc/gamm1; +// const real mach_bc = velocity_magnitude_bc/sound_bc; +// // Eq. 46 +// const real radicant = 1.0+0.5*gamm1*mach_bc*mach_bc; +// const real pressure_bc = total_inlet_pressure * pow(radicant, -gam/gamm1); +// const real temperature_bc = total_inlet_temperature * pow(radicant, -1.0); +// //std::cout << " pressure_bc " << pressure_bc << "pressure_inf" << pressure_inf << std::endl; +// //std::cout << " temperature_bc " << temperature_bc << "temperature_inf" << temperature_inf << std::endl; +// // +// +// const real density_bc = compute_density_from_pressure_temperature(pressure_bc, temperature_bc); +// std::array primitive_boundary_values; +// primitive_boundary_values[0] = density_bc; +// for (int d=0;d primitive_boundary_values; +// primitive_boundary_values[0] = density_bc; +// for (int d=0;d +//dealii::Vector MHD::post_compute_derived_quantities_vector ( +// const dealii::Vector &uh, +// const std::vector > &duh, +// const std::vector > &dduh, +// const dealii::Tensor<1,dim> &normals, +// const dealii::Point &evaluation_points) const +//{ +// std::vector names = post_get_names (); +// dealii::Vector computed_quantities = PhysicsBase::post_compute_derived_quantities_vector ( uh, duh, dduh, normals, evaluation_points); +// unsigned int current_data_index = computed_quantities.size() - 1; +// computed_quantities.grow_or_shrink(names.size()); +// if constexpr (std::is_same::value) { +// +// std::array conservative_soln; +// for (unsigned int s=0; s primitive_soln = convert_conservative_to_primitive(conservative_soln); +// +// // Density +// computed_quantities(++current_data_index) = primitive_soln[0]; +// // Velocities +// for (unsigned int d=0; d +//std::vector MHD +//::post_get_data_component_interpretation () const +//{ +// namespace DCI = dealii::DataComponentInterpretation; +// std::vector interpretation = PhysicsBase::post_get_data_component_interpretation (); // state variables +// interpretation.push_back (DCI::component_is_scalar); // Density +// for (unsigned int d=0; d names = post_get_names(); +// if (names.size() != interpretation.size()) { +// std::cout << "Number of DataComponentInterpretation is not the same as number of names for output file" << std::endl; +// } +// return interpretation; +//} +// +// +//template +//std::vector MHD ::post_get_names () const +//{ +// std::vector names = PhysicsBase::post_get_names (); +// names.push_back ("density"); +// for (unsigned int d=0; d +//dealii::UpdateFlags MHD +//::post_get_needed_update_flags () const +//{ +// //return update_values | update_gradients; +// return dealii::update_values; +//} // Instantiate explicitly -template class MHD < PHILIP_DIM, 2*PHILIP_DIM+2, double >; -template class MHD < PHILIP_DIM, 2*PHILIP_DIM+2, Sacado::Fad::DFad >; +template class MHD < PHILIP_DIM, 8, double >; +template class MHD < PHILIP_DIM, 8, Sacado::Fad::DFad >; } // Physics namespace } // PHiLiP namespace diff --git a/src/physics/mhd.h b/src/physics/mhd.h index d1a47263c..7236b78c0 100644 --- a/src/physics/mhd.h +++ b/src/physics/mhd.h @@ -81,7 +81,7 @@ class MHD : public PhysicsBase : gam(gamma_gas) , gamm1(gam-1.0) { - static_assert(nstate==2*dim+2, "Physics::MHD() should be created with nstate=dim+2"); + static_assert(nstate==8, "Physics::MHD() should be created with nstate=8"); }; /// Destructor @@ -93,6 +93,8 @@ class MHD : public PhysicsBase /// Gamma-1.0 used often const double gamm1; + // double mach_inf_sqr = 1; + std::array manufactured_solution (const dealii::Point &pos) const; /// Convective flux: \f$ \mathbf{F}_{conv} \f$ @@ -207,24 +209,24 @@ class MHD : public PhysicsBase real compute_mean_specific_energy(const std::array &soln_const, const std::array &soln_loop) const; - void boundary_face_values ( - const int /*boundary_type*/, - const dealii::Point &/*pos*/, - const dealii::Tensor<1,dim,real> &/*normal*/, - const std::array &/*soln_int*/, - const std::array,nstate> &/*soln_grad_int*/, - std::array &/*soln_bc*/, - std::array,nstate> &/*soln_grad_bc*/) const; - - virtual dealii::Vector post_compute_derived_quantities_vector ( - const dealii::Vector &uh, - const std::vector > &duh, - const std::vector > &dduh, - const dealii::Tensor<1,dim> &normals, - const dealii::Point &evaluation_points) const; - virtual std::vector post_get_names () const; - virtual std::vector post_get_data_component_interpretation () const; - virtual dealii::UpdateFlags post_get_needed_update_flags () const; +// void boundary_face_values ( +// const int /*boundary_type*/, +// const dealii::Point &/*pos*/, +// const dealii::Tensor<1,dim,real> &/*normal*/, +// const std::array &/*soln_int*/, +// const std::array,nstate> &/*soln_grad_int*/, +// std::array &/*soln_bc*/, +// std::array,nstate> &/*soln_grad_bc*/) const; +// +// virtual dealii::Vector post_compute_derived_quantities_vector ( +// const dealii::Vector &uh, +// const std::vector > &duh, +// const std::vector > &dduh, +// const dealii::Tensor<1,dim> &normals, +// const dealii::Point &evaluation_points) const; +// virtual std::vector post_get_names () const; +// virtual std::vector post_get_data_component_interpretation () const; +// virtual dealii::UpdateFlags post_get_needed_update_flags () const; protected: diff --git a/src/physics/physics.cpp b/src/physics/physics.cpp index a0e8f5b21..7c84c1cf4 100644 --- a/src/physics/physics.cpp +++ b/src/physics/physics.cpp @@ -184,6 +184,8 @@ template class PhysicsBase < PHILIP_DIM, 4, double >; template class PhysicsBase < PHILIP_DIM, 4, Sacado::Fad::DFad >; template class PhysicsBase < PHILIP_DIM, 5, double >; template class PhysicsBase < PHILIP_DIM, 5, Sacado::Fad::DFad >; +template class PhysicsBase < PHILIP_DIM, 8, double >; +template class PhysicsBase < PHILIP_DIM, 8, Sacado::Fad::DFad >; } // Physics namespace } // PHiLiP namespace diff --git a/src/physics/physics_factory.cpp b/src/physics/physics_factory.cpp index b09988841..a0993103a 100644 --- a/src/physics/physics_factory.cpp +++ b/src/physics/physics_factory.cpp @@ -5,7 +5,7 @@ #include "convection_diffusion.h" #include "burgers.h" #include "euler.h" -//#include "mhd.h" +#include "mhd.h" namespace PHiLiP { namespace Physics { @@ -36,9 +36,9 @@ ::create_Physics(const Parameters::AllParameters *const parameters_input) ,parameters_input->euler_param.side_slip_angle); } - } //else if (pde_type == PDE_enum::mhd) { - //if constexpr (nstate == 2*dim + 2) return std::make_shared < MHD > (); - // } + } else if (pde_type == PDE_enum::mhd) { + if constexpr (nstate == 8) return std::make_shared < MHD > (parameters_input->euler_param.gamma_gas); + } std::cout << "Can't create PhysicsBase, invalid PDE type: " << pde_type << std::endl; assert(0==1 && "Can't create PhysicsBase, invalid PDE type"); return nullptr; @@ -54,6 +54,8 @@ template class PhysicsFactory; template class PhysicsFactory >; template class PhysicsFactory; template class PhysicsFactory >; +template class PhysicsFactory; +template class PhysicsFactory >; } // Physics namespace From 8a88c38257c2207079dffcdddf7cdbecc648f63b Mon Sep 17 00:00:00 2001 From: abtin98 Date: Mon, 23 Sep 2019 17:56:18 -0400 Subject: [PATCH 06/10] modified some test parameters --- .../advection_implicit/1d_advection_implicit.prm | 2 ++ .../burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm | 2 ++ .../euler_integration/1d_euler_laxfriedrichs_manufactured.prm | 2 +- .../euler_integration/2d_euler_laxfriedrichs_manufactured.prm | 3 +++ 4 files changed, 8 insertions(+), 1 deletion(-) diff --git a/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit.prm b/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit.prm index af3078458..ef5c35fe3 100644 --- a/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit.prm +++ b/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit.prm @@ -7,6 +7,8 @@ set dimension = 1 # . set pde_type = advection +set use_collocated_nodes = false + set conv_num_flux = lax_friedrichs subsection ODE solver diff --git a/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm b/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm index 1b4b88b45..b4f38831d 100644 --- a/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm +++ b/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm @@ -9,6 +9,8 @@ set pde_type = burgers_inviscid set use_weak_form = true +set use_collocated_nodes = false + subsection ODE solver # Maximum nonlinear solver iterations set nonlinear_max_iterations = 500 diff --git a/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured.prm b/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured.prm index c92dd89ee..af1b0b0e8 100644 --- a/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured.prm +++ b/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured.prm @@ -9,7 +9,7 @@ set conv_num_flux = lax_friedrichs #set use_split_form = true -#set use_collocated_nodes = true +set use_collocated_nodes = false set use_weak_form = true diff --git a/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured.prm b/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured.prm index f8caab53e..967fff158 100644 --- a/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured.prm +++ b/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured.prm @@ -7,6 +7,9 @@ set pde_type = euler set conv_num_flux = lax_friedrichs +set use_weak_form = true + +set use_collocated_nodes = false subsection ODE solver # Maximum nonlinear solver iterations From 423718bc825ce4299dbc6bb1d884790f84d20abb Mon Sep 17 00:00:00 2001 From: abtin98 Date: Mon, 23 Sep 2019 18:09:14 -0400 Subject: [PATCH 07/10] added collocation tests --- .../1d_advection_implicit_collocated.prm | 50 +++++++++++++++ .../advection_implicit/CMakeLists.txt | 8 +++ ...d_burgers_inviscid_implicit_collocated.prm | 45 +++++++++++++ .../burgers_inviscid_implicit/CMakeLists.txt | 8 +++ ...nvection_diffusion_implicit_collocated.prm | 44 +++++++++++++ .../CMakeLists.txt | 7 +++ ..._laxfriedrichs_manufactured_collocated.prm | 57 +++++++++++++++++ ..._laxfriedrichs_manufactured_collocated.prm | 63 +++++++++++++++++++ .../euler_integration/CMakeLists.txt | 16 +++++ 9 files changed, 298 insertions(+) create mode 100644 tests/integration_tests_control_files/advection_implicit/1d_advection_implicit_collocated.prm create mode 100644 tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit_collocated.prm create mode 100644 tests/integration_tests_control_files/convection_diffusion_implicit/2d_convection_diffusion_implicit_collocated.prm create mode 100644 tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured_collocated.prm create mode 100644 tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured_collocated.prm diff --git a/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit_collocated.prm b/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit_collocated.prm new file mode 100644 index 000000000..7d84fc553 --- /dev/null +++ b/tests/integration_tests_control_files/advection_implicit/1d_advection_implicit_collocated.prm @@ -0,0 +1,50 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions +set dimension = 1 + +# The PDE we want to solve. Choices are +# . +set pde_type = advection + +set use_collocated_nodes = true + +set conv_num_flux = lax_friedrichs + +subsection ODE solver + set initial_time_step = 1e10 + # Maximum nonlinear solver iterations + set nonlinear_max_iterations = 500 + + # Nonlinear solver residual tolerance + set nonlinear_steady_residual_tolerance = 1e-13 + + set initial_time_step = 100 + set time_step_factor_residual = 20.0 + set time_step_factor_residual_exp = 3.0 + + # Print every print_iteration_modulo iterations of the nonlinear solver + set print_iteration_modulo = 1 + + # Explicit or implicit solverChoices are . + set ode_solver_type = implicit +end + +subsection manufactured solution convergence study + set use_manufactured_source_term = true + # Last degree used for convergence study + set degree_end = 3 + + # Starting degree for convergence study + set degree_start = 0 + + # Multiplier on grid size. nth-grid will be of size + # (initial_grid^grid_progression)^dim + set grid_progression = 2 + + # Initial grid of size (initial_grid_size)^dim + set initial_grid_size = 2 + + # Number of grids in grid study + set number_of_grids = 4 +end diff --git a/tests/integration_tests_control_files/advection_implicit/CMakeLists.txt b/tests/integration_tests_control_files/advection_implicit/CMakeLists.txt index 0b612e0f4..ce69b9977 100644 --- a/tests/integration_tests_control_files/advection_implicit/CMakeLists.txt +++ b/tests/integration_tests_control_files/advection_implicit/CMakeLists.txt @@ -56,3 +56,11 @@ add_test( COMMAND mpirun -n ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/3d_advection_implicit_strong.prm WORKING_DIRECTORY ${TEST_OUTPUT_DIR} ) + +configure_file(1d_advection_implicit_collocated.prm 1d_advection_implicit_collocated.prm COPYONLY) +add_test( + NAME 1D_ADVECTION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION + COMMAND mpirun -n 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/1d_advection_implicit_collocated.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + diff --git a/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit_collocated.prm b/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit_collocated.prm new file mode 100644 index 000000000..137b9ea6c --- /dev/null +++ b/tests/integration_tests_control_files/burgers_inviscid_implicit/1d_burgers_inviscid_implicit_collocated.prm @@ -0,0 +1,45 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions +set dimension = 1 + +# The PDE we want to solve. Choices are +# . +set pde_type = burgers_inviscid + +set use_weak_form = true + +set use_collocated_nodes = true + +subsection ODE solver + # Maximum nonlinear solver iterations + set nonlinear_max_iterations = 500 + + # Nonlinear solver residual tolerance + set nonlinear_steady_residual_tolerance = 1e-12 + + # Print every print_iteration_modulo iterations of the nonlinear solver + set print_iteration_modulo = 1 + + # Explicit or implicit solverChoices are . + set ode_solver_type = implicit +end + +subsection manufactured solution convergence study + set use_manufactured_source_term = true + # Last degree used for convergence study + set degree_end = 3 + + # Starting degree for convergence study + set degree_start = 1 + + # Multiplier on grid size. nth-grid will be of size + # (initial_grid^grid_progression)^dim + set grid_progression = 1.5 + + # Initial grid of size (initial_grid_size)^dim + set initial_grid_size = 3 + + # Number of grids in grid study + set number_of_grids = 6 +end diff --git a/tests/integration_tests_control_files/burgers_inviscid_implicit/CMakeLists.txt b/tests/integration_tests_control_files/burgers_inviscid_implicit/CMakeLists.txt index 12e49578e..e367a1408 100644 --- a/tests/integration_tests_control_files/burgers_inviscid_implicit/CMakeLists.txt +++ b/tests/integration_tests_control_files/burgers_inviscid_implicit/CMakeLists.txt @@ -7,6 +7,14 @@ add_test( WORKING_DIRECTORY ${TEST_OUTPUT_DIR} ) +configure_file(1d_burgers_inviscid_implicit_collocated.prm 1d_burgers_inviscid_implicit_collocated.prm COPYONLY) +add_test( + NAME 1D_BURGERS_INVISCID_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION + COMMAND mpirun -n 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/1d_burgers_inviscid_implicit_collocated.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + + # configure_file(2d_burgers_inviscid_implicit.prm 2d_burgers_inviscid_implicit.prm COPYONLY) # add_test( # NAME MPI_2D_BURGERS_INVISCID_IMPLICIT_MANUFACTURED_SOLUTION diff --git a/tests/integration_tests_control_files/convection_diffusion_implicit/2d_convection_diffusion_implicit_collocated.prm b/tests/integration_tests_control_files/convection_diffusion_implicit/2d_convection_diffusion_implicit_collocated.prm new file mode 100644 index 000000000..9b3c4647f --- /dev/null +++ b/tests/integration_tests_control_files/convection_diffusion_implicit/2d_convection_diffusion_implicit_collocated.prm @@ -0,0 +1,44 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions +set dimension = 2 + +set use_collocated_nodes = true + +# The PDE we want to solve. Choices are +# . +set pde_type = convection_diffusion + +subsection ODE solver + # Maximum nonlinear solver iterations + set nonlinear_max_iterations = 500 + + # Nonlinear solver residual tolerance + set nonlinear_steady_residual_tolerance = 1e-12 + + # Print every print_iteration_modulo iterations of the nonlinear solver + set print_iteration_modulo = 1 + + # Explicit or implicit solverChoices are . + set ode_solver_type = implicit +end + +subsection manufactured solution convergence study + set use_manufactured_source_term = true + # Last degree used for convergence study + set degree_end = 3 + + # Starting degree for convergence study + set degree_start = 1 + + # Multiplier on grid size. nth-grid will be of size + # (initial_grid^grid_progression)^dim + set grid_progression = 1.5 + + # Initial grid of size (initial_grid_size)^dim + set initial_grid_size = 2 + + # Number of grids in grid study + set number_of_grids = 5 +end + diff --git a/tests/integration_tests_control_files/convection_diffusion_implicit/CMakeLists.txt b/tests/integration_tests_control_files/convection_diffusion_implicit/CMakeLists.txt index b14c7d99c..3ced36f8a 100644 --- a/tests/integration_tests_control_files/convection_diffusion_implicit/CMakeLists.txt +++ b/tests/integration_tests_control_files/convection_diffusion_implicit/CMakeLists.txt @@ -14,6 +14,13 @@ add_test( WORKING_DIRECTORY ${TEST_OUTPUT_DIR} ) +configure_file(2d_convection_diffusion_implicit_collocated.prm 2d_convection_diffusion_implicit_collocated.prm COPYONLY) +add_test( + NAME 2D_CONVECTION_DIFFUSION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION + COMMAND mpirun -n 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/2d_convection_diffusion_implicit_collocated.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + configure_file(3d_convection_diffusion_implicit.prm 3d_convection_diffusion_implicit.prm COPYONLY) add_test( NAME MPI_3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION diff --git a/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured_collocated.prm b/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured_collocated.prm new file mode 100644 index 000000000..c5df825dd --- /dev/null +++ b/tests/integration_tests_control_files/euler_integration/1d_euler_laxfriedrichs_manufactured_collocated.prm @@ -0,0 +1,57 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions +set dimension = 1 + +set pde_type = euler + +set conv_num_flux = lax_friedrichs + +#set use_split_form = true + +set use_collocated_nodes = true + +set use_weak_form = true + +subsection ODE solver + + set ode_output = verbose + + set initial_time_step = 1000 + set time_step_factor_residual = 10 + set time_step_factor_residual_exp = 2 + + # Maximum nonlinear solver iterations + set nonlinear_max_iterations = 500000 + + # Nonlinear solver residual tolerance + set nonlinear_steady_residual_tolerance = 1e-12 + + # Print every print_iteration_modulo iterations of the nonlinear solver + set print_iteration_modulo = 1 + + # Explicit or implicit solverChoices are . + set ode_solver_type = implicit +end + +subsection manufactured solution convergence study + set use_manufactured_source_term = true + # Last degree used for convergence study + set degree_end = 3 + + # Starting degree for convergence study + set degree_start = 0 + + # Multiplier on grid size. nth-grid will be of size + # (initial_grid^grid_progression)^dim + set grid_progression = 2 + + set grid_progression_add = 5 + # Initial grid of size (initial_grid_size)^dim + set initial_grid_size = 10 + + # Number of grids in grid study + set number_of_grids = 4 + + set slope_deficit_tolerance = 0.2 +end diff --git a/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured_collocated.prm b/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured_collocated.prm new file mode 100644 index 000000000..b4f612b6b --- /dev/null +++ b/tests/integration_tests_control_files/euler_integration/2d_euler_laxfriedrichs_manufactured_collocated.prm @@ -0,0 +1,63 @@ +# Listing of Parameters +# --------------------- +# Number of dimensions +set dimension = 2 + +set pde_type = euler + +set conv_num_flux = lax_friedrichs + +set use_weak_form = true + +set use_collocated_nodes = true + +subsection ODE solver + # Maximum nonlinear solver iterations + set nonlinear_max_iterations = 100 + + # Nonlinear solver residual tolerance + set nonlinear_steady_residual_tolerance = 1e-9 + + set initial_time_step = 1000 + set time_step_factor_residual = 20.0 + set time_step_factor_residual_exp = 2.0 + + # Print every print_iteration_modulo iterations of the nonlinear solver + set print_iteration_modulo = 1 + + # Explicit or implicit solverChoices are . + set ode_solver_type = implicit +end + +subsection linear solver + subsection gmres options + set max_iterations = 200 + set linear_residual_tolerance = 1e-4 + set restart_number = 60 + end +end + +subsection manufactured solution convergence study + set use_manufactured_source_term = true + # Last degree used for convergence study + set degree_end = 3 + + # Starting degree for convergence study + set degree_start = 0 + + set grid_progression = 1.0 + + set grid_progression_add = 5 + + # Initial grid of size (initial_grid_size)^dim + set initial_grid_size = 5 + + # Number of grids in grid study + set number_of_grids = 4 + + # WARNING + # If we want actual optimal orders with a tigher tolerance + # we need to increase the grid sizes by a significant amount + set slope_deficit_tolerance = 0.1 +end + diff --git a/tests/integration_tests_control_files/euler_integration/CMakeLists.txt b/tests/integration_tests_control_files/euler_integration/CMakeLists.txt index 3700152ff..e1a123359 100644 --- a/tests/integration_tests_control_files/euler_integration/CMakeLists.txt +++ b/tests/integration_tests_control_files/euler_integration/CMakeLists.txt @@ -10,6 +10,14 @@ add_test( WORKING_DIRECTORY ${TEST_OUTPUT_DIR} ) +configure_file(1d_euler_laxfriedrichs_manufactured_collocated.prm 1d_euler_laxfriedrichs_manufactured_collocated.prm COPYONLY) +add_test( + NAME 1D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION + COMMAND mpirun -np 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/1d_euler_laxfriedrichs_manufactured_collocated.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + + configure_file(2d_euler_laxfriedrichs_manufactured.prm 2d_euler_laxfriedrichs_manufactured.prm COPYONLY) add_test( NAME 2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION @@ -17,6 +25,14 @@ add_test( WORKING_DIRECTORY ${TEST_OUTPUT_DIR} ) +configure_file(2d_euler_laxfriedrichs_manufactured_collocated.prm 2d_euler_laxfriedrichs_manufactured_collocated.prm COPYONLY) +add_test( + NAME 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION + COMMAND mpirun -np 1 ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/2d_euler_laxfriedrichs_manufactured_collocated.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + + configure_file(2d_euler_laxfriedrichs_manufactured.prm 2d_euler_laxfriedrichs_manufactured.prm COPYONLY) add_test( NAME MPI_2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION From d8e784b70d9c799a10ea04d9cfa2ead4ffd964cb Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Wed, 2 Oct 2019 00:05:58 -0400 Subject: [PATCH 08/10] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 409ac0a99..15b24706a 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ - Supported Partial Differential Equations: Linear advection, diffusion, convection-diffusion, Burgers, Euler, TODO: Navier-Stokes. - Supported convective numerical fluxes: Lax-Friedrichs, Roe (Harten's entropy fix) for Euler, InProgress: Split-Form - Supported diffusive numerical fluxes: Symmetric Interior Penalty -- Supported elements: LINEs, QUADs, HEXs +- Supported elements: LINEs, QUADs, HEXs since it uses deal.II - Supported refinements: h (size) or p (order). ## Building/Running the Code From 3f60a58006578f20fd1508100fa07315ac442968 Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Wed, 2 Oct 2019 00:06:48 -0400 Subject: [PATCH 09/10] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 15b24706a..8db1c1cb7 100644 --- a/README.md +++ b/README.md @@ -99,7 +99,7 @@ GDB$ quit If the error only occurs when using parallelism, you can use the following example command ```sh -mpirun -np 2 xterm -hold -e gdb -ex 'break MPI_Abort' -ex run --args /home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/bin/PHiLiP_2D "-i" "/home/ddong/Codes/PHiLiP_temp/PHiLiP/build_de bug/tests/advection_implicit/2d_advection_implicit_strong.prm" +mpirun -np 2 xterm -hold -e gdb -ex 'break MPI_Abort' -ex run --args /home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/bin/PHiLiP_2D "-i" "/home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/tests/advection_implicit/2d_advection_implicit_strong.prm" ``` This launches 2 xterm processes, each of which will launch gdb processes that will run the code and will have a breakpoint when MPI_Abort is encountered. From acccccc4608dff9d174a5fba68e6684672a4208d Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Wed, 2 Oct 2019 00:09:38 -0400 Subject: [PATCH 10/10] Update README.md --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 8db1c1cb7..482065553 100644 --- a/README.md +++ b/README.md @@ -70,6 +70,7 @@ ROOT$ ctest -R (Run tests matching regular expression) ROOT$ ctest -E (Exclude tests matching regular expression) ROOT$ ctest -V (Enable verbose output from tests) ``` +Note that running `ctest` in `Debug` will take forever since some integration tests fully solve nonlinear problems with multiple orders and multiple meshes. It is suggested to perform `ctest` in `Release` mode, and only use `Debug` mode for debugging purposes. ## Debugging