From 4dd6d62e371e262c254dafb377fa536a6d26a413 Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Mon, 23 Sep 2019 18:31:33 -0400 Subject: [PATCH] Add basic HighOrderGrid class and test. --- src/CMakeLists.txt | 14 +- src/dg/CMakeLists.txt | 1 + src/dg/dg.cpp | 317 ++++++++++++++---- src/dg/dg.h | 175 +++++++--- src/dg/high_order_grid.cpp | 47 +++ src/dg/high_order_grid.h | 71 ++++ src/dg/strong_dg.cpp | 21 +- src/dg/weak_dg.cpp | 17 +- src/testing/burgers_stability.cpp | 3 +- src/testing/euler_cylinder.cpp | 3 +- src/testing/euler_entropy_waves.cpp | 3 +- src/testing/euler_gaussian_bump.cpp | 3 +- src/testing/euler_vortex.cpp | 3 +- src/testing/grid_study.cpp | 6 +- .../advection_explicit_periodic.cpp | 3 +- ...ler_split_inviscid_taylor_green_vortex.cpp | 3 +- tests/unit_tests/CMakeLists.txt | 1 + tests/unit_tests/grid/CMakeLists.txt | 36 ++ .../unit_tests/grid/high_order_grid_test.cpp | 227 +++++++++++++ .../regression/jacobian_matrix_regression.cpp | 3 +- 20 files changed, 812 insertions(+), 145 deletions(-) create mode 100644 src/dg/high_order_grid.cpp create mode 100644 src/dg/high_order_grid.h create mode 100644 tests/unit_tests/grid/CMakeLists.txt create mode 100644 tests/unit_tests/grid/high_order_grid_test.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b2b7e1e28..b3623d695 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -76,16 +76,24 @@ set(SMALL_TEST testrun) set(TARGET ${SMALL_TEST}) set(SMALL_TEST_SRC temporary_test.cpp) - add_executable(${SMALL_TEST} ${SMALL_TEST_SRC}) target_link_libraries(${SMALL_TEST} mpi ) - DEAL_II_SETUP_TARGET(${SMALL_TEST}) message("Adding executable " ${SMALL_TEST} " with files " ${SMALL_TEST_SRC} "\n") - unset(SMALL_TEST) unset(TARGET ${SMALL_TEST}) +#set(SMALL_TEST instantiate_vector_ad) +#set(TARGET ${SMALL_TEST}) +#set(SMALL_TEST_SRC +# instantiate_vector_ad.cpp) +#add_executable(${SMALL_TEST} ${SMALL_TEST_SRC}) +#target_link_libraries(${SMALL_TEST} mpi ) +#DEAL_II_SETUP_TARGET(${SMALL_TEST}) +#message("Adding executable " ${SMALL_TEST} " with files " ${SMALL_TEST_SRC} "\n") +#unset(SMALL_TEST) +#unset(TARGET ${SMALL_TEST}) + set(SMALL_TEST manifoldrun) set(TARGET ${SMALL_TEST}) set(SMALL_TEST_SRC diff --git a/src/dg/CMakeLists.txt b/src/dg/CMakeLists.txt index 715fe3b56..0ceb8abf2 100644 --- a/src/dg/CMakeLists.txt +++ b/src/dg/CMakeLists.txt @@ -1,4 +1,5 @@ set(DG_SOURCE + high_order_grid.cpp dg.cpp weak_dg.cpp strong_dg.cpp diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp index f4ad0077d..3e24192c6 100644 --- a/src/dg/dg.cpp +++ b/src/dg/dg.cpp @@ -25,17 +25,156 @@ //#include // Might need mapping_q #include // Might need mapping_q #include +#include // Finally, we take our exact solution from the library as well as volume_quadrature // and additional tools. #include #include +#include +#include #include "dg.h" #include "post_processor/physics_post_processor.h" +//template class dealii::MappingFEField, dealii::hp::DoFHandler >; namespace PHiLiP { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + template using Triangulation = dealii::Triangulation; +#else + template using Triangulation = dealii::parallel::distributed::Triangulation; +#endif + +//template +// > +//void +//get_position_vector2(const dealii::hp::DoFHandler &dh, VectorType &vector) +//{ +// AssertDimension(vector.size(), dh.n_dofs()); +// const dealii::hp::FECollection &fe_coll = dh.get_fe_collection(); +// +// dealii::hp::MappingCollection mapping_coll; +// dealii::hp::QCollection quad_coll; +// for (int i=0; i finite_element = fe_coll[i]; +// +// Assert(finite_element.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); +// assert(finite_element.has_support_points()); +// AssertDimension(finite_element.dofs_per_cell, finite_element.get_unit_support_points().size()); +// +// const dealii::Quadrature quadrature(finite_element.get_unit_support_points()); +// const dealii::MappingQ mapping(finite_element.degree); +// +// mapping_coll.push_back(mapping); +// quad_coll.push_back(quadrature); +// } +// +// dealii::hp::FEValues fe_values_coll (mapping_coll, fe_coll, quadrature_coll, update_quadrature_points); +// +// // Construct default fe_mask; +// const ComponentMask fe_mask(//mask.size() ? mask : +// ComponentMask(fe_coll.get_nonzero_components(0).size(), true)); +// +// AssertDimension(fe_mask.size(), fe_coll.get_nonzero_components(0).size()); +// +// std::vector fe_to_real(fe_mask.size(), +// numbers::invalid_unsigned_int); +// unsigned int size = 0; +// for (unsigned int i = 0; i < fe_mask.size(); ++i) { +// if (fe_mask[i]) fe_to_real[i] = size++; +// } +// Assert(size == spacedim, +// ExcMessage("The Component Mask you provided is invalid. It has to select exactly spacedim entries.")); +// +// std::vector dof_indices(fe_coll[fe_coll.size()-1].dofs_per_cell); +// for (const auto &cell : dh.active_cell_iterators()) { +// if (cell->is_locally_owned()) { +// const unsigned int fe_index = cell->active_fe_index(); +// fe_values_coll.reinit (cell, fe_index, fe_index, fe_index); +// dof_indices.resize(fe_coll[fe_index].dofs_per_cell); +// cell->get_dof_indices(dof_indices); +// const std::vector> &points = fe_v.get_quadrature_points(); +// for (unsigned int q = 0; q < points.size(); ++q) { +// +// const unsigned int comp = fe_coll.system_to_component_index(q).first; +// if (fe_mask[comp]) { +// dealii::internal::ElementAccess:: +// set(points[q][fe_to_real[comp]], dof_indices[q], vector); +// } +// } +// } +// } +//} + +//template +// > +//void +//get_position_vector2(const dealii::hp::DoFHandler &dof_handler, VectorType &vector) +//{ +// AssertDimension(vector.size(), dof_handler.n_dofs()); +// const dealii::hp::FECollection &fe_coll = dof_handler.get_fe_collection(); +// +// dealii::hp::MappingCollection mapping_coll; +// dealii::hp::QCollection quad_coll; +// for (unsigned int i=0; i &finite_element = fe_coll[i]; +// +// Assert(finite_element.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); +// assert(finite_element.has_support_points()); +// AssertDimension(finite_element.dofs_per_cell, finite_element.get_unit_support_points().size()); +// +// const dealii::Quadrature quadrature(finite_element.get_unit_support_points()); +// const dealii::MappingQ mapping(finite_element.degree); +// +// mapping_coll.push_back(mapping); +// quad_coll.push_back(quadrature); +// } +// +// dealii::hp::FEValues fe_values_coll (mapping_coll, fe_coll, quad_coll, dealii::update_quadrature_points); +// +// // Construct default fe_mask; +// const dealii::ComponentMask fe_mask(//mask.size() ? mask : +// dealii::ComponentMask(fe_coll[0].get_nonzero_components(0).size(), true)); +// +// AssertDimension(fe_mask.size(), fe_coll[0].get_nonzero_components(0).size()); +// +// std::vector fe_to_real(fe_mask.size(), +// dealii::numbers::invalid_unsigned_int); +// unsigned int size = 0; +// for (unsigned int i = 0; i < fe_mask.size(); ++i) { +// if (fe_mask[i]) fe_to_real[i] = size++; +// } +// Assert(size == spacedim, +// ExcMessage("The Component Mask you provided is invalid. It has to select exactly spacedim entries.")); +// +// std::vector dof_indices(fe_coll[fe_coll.size()-1].dofs_per_cell); +// for (const auto &cell : dof_handler.active_cell_iterators()) { +// if (cell->is_locally_owned()) { +// const unsigned int fe_index = cell->active_fe_index(); +// fe_values_coll.reinit (cell, fe_index, fe_index, fe_index); +// const dealii::FEValues &fe_values_volume = fe_values_coll.get_present_fe_values(); +// +// dof_indices.resize(fe_coll[fe_index].dofs_per_cell); +// cell->get_dof_indices(dof_indices); +// +// const std::vector> &points = fe_values_volume.get_quadrature_points(); +// for (unsigned int q = 0; q < points.size(); ++q) { +// +// const unsigned int comp = fe_values_volume.get_fe().system_to_component_index(q).first; +// if (fe_mask[comp]) { +// dealii::internal::ElementAccess:: +// set(points[q][fe_to_real[comp]], dof_indices[q], vector); +// } +// } +// } +// } +//} + // DGFactory *********************************************************************** template @@ -43,38 +182,39 @@ std::shared_ptr< DGBase > DGFactory ::create_discontinuous_galerkin( const Parameters::AllParameters *const parameters_input, - const unsigned int degree) + const unsigned int degree, + Triangulation *const triangulation_input) { using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation; PDE_enum pde_type = parameters_input->pde_type; if (parameters_input->use_weak_form) { if (pde_type == PDE_enum::advection) { - return std::make_shared< DGWeak >(parameters_input, degree); + return std::make_shared< DGWeak >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::advection_vector) { - return std::make_shared< DGWeak >(parameters_input, degree); + return std::make_shared< DGWeak >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::diffusion) { - return std::make_shared< DGWeak >(parameters_input, degree); + return std::make_shared< DGWeak >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::convection_diffusion) { - return std::make_shared< DGWeak >(parameters_input, degree); + return std::make_shared< DGWeak >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::burgers_inviscid) { - return std::make_shared< DGWeak >(parameters_input, degree); + return std::make_shared< DGWeak >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::euler) { - return std::make_shared< DGWeak >(parameters_input, degree); + return std::make_shared< DGWeak >(parameters_input, degree, triangulation_input); } } else { if (pde_type == PDE_enum::advection) { - return std::make_shared< DGStrong >(parameters_input, degree); + return std::make_shared< DGStrong >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::advection_vector) { - return std::make_shared< DGStrong >(parameters_input, degree); + return std::make_shared< DGStrong >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::diffusion) { - return std::make_shared< DGStrong >(parameters_input, degree); + return std::make_shared< DGStrong >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::convection_diffusion) { - return std::make_shared< DGStrong >(parameters_input, degree); + return std::make_shared< DGStrong >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::burgers_inviscid) { - return std::make_shared< DGStrong >(parameters_input, degree); + return std::make_shared< DGStrong >(parameters_input, degree, triangulation_input); } else if (pde_type == PDE_enum::euler) { - return std::make_shared< DGStrong >(parameters_input, degree); + return std::make_shared< DGStrong >(parameters_input, degree, triangulation_input); } } std::cout << "Can't create DGBase in create_discontinuous_galerkin(). Invalid PDE type: " << pde_type << std::endl; @@ -86,8 +226,9 @@ template DGBase::DGBase( // @suppress("Class members should be properly initialized") const int nstate_input, const Parameters::AllParameters *const parameters_input, - const unsigned int max_degree_input) - : DGBase(nstate_input, parameters_input, max_degree_input, this->create_collection_tuple(max_degree_input, nstate_input, parameters_input)) + const unsigned int max_degree_input, + Triangulation *const triangulation_input) + : DGBase(nstate_input, parameters_input, max_degree_input, triangulation_input, this->create_collection_tuple(max_degree_input, nstate_input, parameters_input)) { } template @@ -95,35 +236,53 @@ DGBase::DGBase( // @suppress("Class members should be properly initial const int nstate_input, const Parameters::AllParameters *const parameters_input, const unsigned int max_degree_input, - std::tuple< dealii::hp::MappingCollection, dealii::hp::FECollection, - dealii::hp::QCollection, dealii::hp::QCollection, dealii::hp::QCollection<1>, - dealii::hp::FECollection > collection_tuple) + Triangulation *const triangulation_input, + const MassiveCollectionTuple collection_tuple) : all_parameters(parameters_input) , nstate(nstate_input) , max_degree(max_degree_input) - , mapping_collection(std::get<0>(collection_tuple)) - , fe_collection(std::get<1>(collection_tuple)) + , triangulation(triangulation_input) + //, mapping_collection(std::get<0>(collection_tuple)) + , fe_collection(std::get<0>(collection_tuple)) + , fe_grid(std::get<1>(collection_tuple),dim) , volume_quadrature_collection(std::get<2>(collection_tuple)) , face_quadrature_collection(std::get<3>(collection_tuple)) , oned_quadrature_collection(std::get<4>(collection_tuple)) - , fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags) - , fe_values_collection_face_int (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags) - , fe_values_collection_face_ext (mapping_collection, fe_collection, face_quadrature_collection, this->neighbor_face_update_flags) - , fe_values_collection_subface (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags) , fe_collection_lagrange(std::get<5>(collection_tuple)) - , fe_values_collection_volume_lagrange (mapping_collection, fe_collection_lagrange, volume_quadrature_collection, this->volume_update_flags) + , dof_handler(*triangulation) + , dof_handler_grid(*triangulation) + , high_order_grid_mapping(dof_handler_grid, high_order_grid) , mpi_communicator(MPI_COMM_WORLD) , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) -{} +{ -template -std::tuple< dealii::hp::MappingCollection, dealii::hp::FECollection, - dealii::hp::QCollection, dealii::hp::QCollection, dealii::hp::QCollection<1>, - dealii::hp::FECollection > + dof_handler.initialize(*triangulation, fe_collection); + dof_handler_grid.initialize(*triangulation, fe_grid); + int minimum_degree = (all_parameters->use_collocated_nodes==true) ? 1 : 0; + for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) { + mapping_collection.push_back(high_order_grid_mapping); // Can't just return FESystem. For some reason the copy constructor is deleted + } + + set_all_cells_fe_degree (fe_collection.size()-1); // Always sets it to the maximum degree + +} + +template +std::tuple< + //dealii::hp::MappingCollection, // Mapping + dealii::hp::FECollection, // Solution FE + //dealii::hp::FECollection, // Grid FE + dealii::FE_Q, // Grid FE + dealii::hp::QCollection, // Volume quadrature + dealii::hp::QCollection, // Face quadrature + dealii::hp::QCollection<1>, // 1D quadrature for strong form + dealii::hp::FECollection > // Lagrange polynomials for strong form DGBase::create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const { - dealii::hp::MappingCollection mapping_coll; + //dealii::hp::MappingCollection mapping_coll; dealii::hp::FECollection fe_coll; + //dealii::hp::FECollection fe_coll_grid; + const dealii::FE_Q feq_grid(max_degree+1); // The grid must be at least p1. A p0 solution required a p1 grid. dealii::hp::QCollection volume_quad_coll; dealii::hp::QCollection face_quad_coll; dealii::hp::QCollection<1> oned_quad_coll; @@ -133,14 +292,19 @@ DGBase::create_collection_tuple(const unsigned int max_degree, const i for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) { //const dealii::MappingQ mapping(degree, true); //const dealii::MappingQ mapping(degree+1, true); - const dealii::MappingManifold mapping; - mapping_coll.push_back(mapping); + //const dealii::MappingManifold mapping; + //mapping_coll.push_back(mapping); + //mapping_coll.push_back(high_order_grid_mapping); // Can't just return FESystem. For some reason the copy constructor is deleted + // Solution FECollection const dealii::FE_DGQ fe_dg(degree); const dealii::FESystem fe_system(fe_dg, nstate); fe_coll.push_back (fe_system); - // + //// Grid FECollection uses Lagrange polynomials + //const dealii::FE_Q fe_q(degree+1); // The grid must be at least p1. A p0 solution required a p1 grid. + //const dealii::FESystem fe_system_q(fe_q, dim); + //fe_coll_grid.push_back (fe_system_q); dealii::Quadrature<1> oned_quad(degree+1); dealii::Quadrature volume_quad(degree+1); @@ -185,7 +349,7 @@ DGBase::create_collection_tuple(const unsigned int max_degree, const i dealii::FE_DGQArbitraryNodes lagrange_poly(oned_quad); fe_coll_lagr.push_back (lagrange_poly); } - return std::make_tuple(mapping_coll, fe_coll, volume_quad_coll, face_quad_coll, oned_quad_coll, fe_coll_lagr); + return std::make_tuple(fe_coll, feq_grid, volume_quad_coll, face_quad_coll, oned_quad_coll, fe_coll_lagr); } @@ -197,6 +361,11 @@ void DGBase::set_all_cells_fe_degree ( const unsigned int degree ) { if (cell->is_locally_owned()) cell->set_future_fe_index (degree); } + + //for (auto cell = dof_handler_grid.begin_active(); cell != dof_handler_grid.end(); ++cell) + //{ + // if (cell->is_locally_owned()) cell->set_future_fe_index (degree); + //} triangulation->execute_coarsening_and_refinement(); } @@ -207,24 +376,24 @@ template DGBase::~DGBase () { dof_handler.clear (); + dof_handler_grid.clear (); } -#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D - template using Triangulation = dealii::Triangulation; -#else - template using Triangulation = dealii::parallel::distributed::Triangulation; -#endif -template -void DGBase::set_triangulation(Triangulation *triangulation_input) -{ - dof_handler.clear(); - triangulation = triangulation_input; - dof_handler.initialize(*triangulation, fe_collection); - - set_all_cells_fe_degree (fe_collection.size()-1); // Always sets it to the maximum degree - - //set_all_cells_fe_degree ( max_degree-min_degree); -} +// template +// void DGBase::set_triangulation(Triangulation *triangulation_input) +// { +// // dof_handler.clear(); +// // dof_handler_grid.clear(); +// +// // triangulation = triangulation_input; +// +// // dof_handler.initialize(*triangulation, fe_collection); +// // dof_handler_grid.initialize(*triangulation, fe_grid); +// +// //set_all_cells_fe_degree (fe_collection.size()-1); // Always sets it to the maximum degree +// +// //set_all_cells_fe_degree ( max_degree-min_degree); +// } template @@ -239,6 +408,13 @@ void DGBase::assemble_residual (const bool compute_dRdW) std::vector current_dofs_indices(max_dofs_per_cell); std::vector neighbor_dofs_indices(max_dofs_per_cell); + dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. + dealii::hp::FEFaceValues fe_values_collection_face_int (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of interior face. + dealii::hp::FEFaceValues fe_values_collection_face_ext (mapping_collection, fe_collection, face_quadrature_collection, this->neighbor_face_update_flags); ///< FEValues of exterior face. + dealii::hp::FESubfaceValues fe_values_collection_subface (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of subface. + + dealii::hp::FEValues fe_values_collection_volume_lagrange (mapping_collection, fe_collection_lagrange, volume_quadrature_collection, this->volume_update_flags); + unsigned int n_cell_visited = 0; unsigned int n_face_visited = 0; @@ -264,13 +440,12 @@ void DGBase::assemble_residual (const bool compute_dRdW) fe_values_collection_volume.reinit (current_cell, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); - if (!(all_parameters->use_weak_form)) { - fe_values_collection_volume_lagrange.reinit (current_cell, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); - } + if (!(all_parameters->use_weak_form)) fe_values_collection_volume_lagrange.reinit (current_cell, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + const dealii::FEValues &fe_values_lagrange = fe_values_collection_volume_lagrange.get_present_fe_values(); if ( compute_dRdW ) { - assemble_volume_terms_implicit (fe_values_volume, current_dofs_indices, current_cell_rhs); + assemble_volume_terms_implicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); } else { - assemble_volume_terms_explicit (fe_values_volume, current_dofs_indices, current_cell_rhs); + assemble_volume_terms_explicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); } for (unsigned int iface=0; iface < dealii::GeometryInfo::faces_per_cell; ++iface) { @@ -692,15 +867,16 @@ void DGBase::output_results_vtk (const unsigned int cycle)// const data_out.add_data_vector (right_hand_side, residual_names, dealii::DataOut_DoFData,dim>::DataVectorType::type_dof_data); + const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); data_out.build_patches (mapping_collection[mapping_collection.size()-1]); std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; filename += dealii::Utilities::int_to_string(cycle, 4) + "."; - filename += dealii::Utilities::int_to_string(triangulation->locally_owned_subdomain(), 4); + filename += dealii::Utilities::int_to_string(iproc, 4); filename += ".vtu"; std::ofstream output(filename); data_out.write_vtu(output); - if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) == 0) { + if (iproc == 0) { std::vector filenames; for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) { std::string fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; @@ -717,7 +893,6 @@ void DGBase::output_results_vtk (const unsigned int cycle)// const } - template void DGBase::allocate_system () { @@ -726,6 +901,29 @@ void DGBase::allocate_system () // system matrices and vectors. dof_handler.distribute_dofs(fe_collection); + dof_handler_grid.distribute_dofs(fe_grid); + + // High order grid initialization + locally_owned_dofs_grid = dof_handler_grid.locally_owned_dofs(); + dealii::DoFTools::extract_locally_relevant_dofs(dof_handler_grid, ghost_dofs_grid); + locally_relevant_dofs_grid = ghost_dofs_grid; + ghost_dofs_grid.subtract_set(locally_owned_dofs_grid); + high_order_grid.reinit(locally_owned_dofs_grid, ghost_dofs_grid, mpi_communicator); + + //dealii::VectorTools::get_position_vector2(dof_handler_grid, high_order_grid); + //get_position_vector2(dof_handler_grid, high_order_grid); + dealii::MappingFEField, dealii::DoFHandler > map(dof_handler_grid, high_order_grid); + //int minimum_degree = (all_parameters->use_collocated_nodes==true) ? 1 : 0; + ////int current_fe_index = 0; + //for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) { + // mapping_collection.push_back(map); + // //if(current_fe_index <= mapping_collection.size()) { + // // mapping_collection.push_back(map); + // //} else { + // // mapping_collection[current_fe_index] = map; + // //} + // //current_fe_index++; + //} // Solution and RHS locally_owned_dofs = dof_handler.locally_owned_dofs(); @@ -796,6 +994,7 @@ void DGBase::evaluate_mass_matrices (bool do_inverse_mass_matrix) //} //pcout << "AFter reinit" << std::endl; + dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) { if (!cell->is_locally_owned()) continue; diff --git a/src/dg/dg.h b/src/dg/dg.h index 922f7ea15..5dfc78a66 100644 --- a/src/dg/dg.h +++ b/src/dg/dg.h @@ -9,6 +9,7 @@ #include #include #include +#include #include @@ -28,8 +29,22 @@ #include "numerical_flux/numerical_flux.h" #include "parameters/all_parameters.h" +// Template specialization of MappingFEField +//extern template class dealii::MappingFEField, dealii::hp::DoFHandler >; namespace PHiLiP { +//#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +// using Triangulation = dealii::Triangulation; +//#else +// using Triangulation = dealii::parallel::distributed::Triangulation; +//#endif +//namespace PHiLiP { +//#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +// template using Triangulation = dealii::Triangulation; +//#else +// template using Triangulation = dealii::parallel::distributed::Triangulation; +//#endif + /// DGBase is independent of the number of state variables. /** This base class allows the use of arrays to efficiently allocate the data structures * through std::array in the derived class DG. @@ -50,6 +65,11 @@ namespace PHiLiP { template class DGBase { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + using Triangulation = dealii::Triangulation; +#else + using Triangulation = dealii::parallel::distributed::Triangulation; +#endif public: const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters @@ -70,7 +90,23 @@ class DGBase * * Passes create_collection_tuple() to the delegated constructor. */ - DGBase(const int nstate_input, const Parameters::AllParameters *const parameters_input, const unsigned int max_degree); + DGBase(const int nstate_input, + const Parameters::AllParameters *const parameters_input, + const unsigned int max_degree, + Triangulation *const triangulation_input); + + + /// Makes for cleaner doxygen documentation + using MassiveCollectionTuple = std::tuple< + //dealii::hp::MappingCollection, // Mapping + dealii::hp::FECollection, // Solution FE + //dealii::hp::FECollection, // Grid FE + //dealii::FESystem, // Grid FE // Can't just return FESystem. For some reason the copy constructor is deleted + dealii::FE_Q, // Grid FE + dealii::hp::QCollection, // Volume quadrature + dealii::hp::QCollection, // Face quadrature + dealii::hp::QCollection<1>, // 1D quadrature for strong form + dealii::hp::FECollection >; // Lagrange polynomials for strong form /// Delegated constructor that initializes collections. /** Since a function is used to generate multiple different objects, a delegated @@ -80,20 +116,20 @@ class DGBase DGBase( const int nstate_input, const Parameters::AllParameters *const parameters_input, const unsigned int max_degree_input, - const std::tuple< dealii::hp::MappingCollection, dealii::hp::FECollection, - dealii::hp::QCollection, dealii::hp::QCollection, dealii::hp::QCollection<1>, - dealii::hp::FECollection > collection_tuple); + Triangulation *const triangulation_input, + const MassiveCollectionTuple collection_tuple); virtual ~DGBase(); ///< Destructor. -#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D - using Triangulation = dealii::Triangulation; -#else - using Triangulation = dealii::parallel::distributed::Triangulation; -#endif - Triangulation *triangulation; ///< Mesh + Triangulation *const triangulation; ///< Mesh /// Sets the triangulation for 2D and 3D. Should be done before allocate system - void set_triangulation(Triangulation *triangulation_input); + //void set_triangulation(Triangulation *triangulation_input); + + /// Refers to a collection Mappings, which represents the high-order grid. + /** Since we are interested in performing mesh movement for optimization purposes, + * this is not a constant member variables. + */ + dealii::hp::MappingCollection mapping_collection; void set_all_cells_fe_degree ( const unsigned int degree ); @@ -124,16 +160,6 @@ class DGBase unsigned int n_dofs() const; ///< Number of degrees of freedom - /// Degrees of freedom handler - /* Allows us to iterate over the finite elements' degrees of freedom. - * Note that since we are not using FESystem, we need to multiply - * the index by a factor of "nstate" - * - * Must be defined after fe_dg since it is a subscriptor of fe_dg. - * Destructor are called in reverse order in which they appear in class definition. - */ - dealii::hp::DoFHandler dof_handler; - /// Sparsity pattern used on the system_matrix /** Not sure we need to store it. */ dealii::SparsityPattern sparsity_pattern; @@ -183,12 +209,19 @@ class DGBase dealii::IndexSet locally_owned_dofs; ///< Locally own degrees of freedom dealii::IndexSet ghost_dofs; ///< Locally relevant ghost degrees of freedom dealii::IndexSet locally_relevant_dofs; ///< Union of locally owned degrees of freedom and relevant ghost degrees of freedom + + dealii::IndexSet locally_owned_dofs_grid; ///< Locally own degrees of freedom for the grid + dealii::IndexSet ghost_dofs_grid; ///< Locally relevant ghost degrees of freedom for the grid + dealii::IndexSet locally_relevant_dofs_grid; ///< Union of locally owned degrees of freedom and relevant ghost degrees of freedom for the grid /// Current modal coefficients of the solution /** Note that the current processor has read-access to all locally_relevant_dofs * and has write-access to all locally_owned_dofs */ dealii::LinearAlgebra::distributed::Vector solution; + /// Current nodal coefficients of the high-order grid. + dealii::LinearAlgebra::distributed::Vector high_order_grid; + void initialize_manufactured_solution (); ///< Virtual function defined in DG void output_results_vtk (const unsigned int ith_grid); ///< Output solution @@ -228,15 +261,18 @@ class DGBase //void assemble_residual_dRdW (); void assemble_residual (const bool compute_dRdW=false); - /// Mapping is currently MappingQ. - /** Refer to deal.II documentation for the various mapping types */ - //const dealii::MappingQ mapping; - const dealii::hp::MappingCollection mapping_collection; - - /// Finite Element Collection for p-finite-element + /// Finite Element Collection for p-finite-element to represent the solution /** This is a collection of FESystems */ const dealii::hp::FECollection fe_collection; + /// Finite Element Collection to represent the high-order grid + /** This is a collection of FESystems. + * Unfortunately, deal.II doesn't have a working hp Mapping FE field. + * Therefore, every grid/cell will use the maximal polynomial mapping regardless of the solution order. + */ + //const dealii::hp::FECollection fe_collection_grid; + const dealii::FESystem fe_grid; + dealii::hp::QCollection volume_quadrature_collection; dealii::hp::QCollection face_quadrature_collection; dealii::hp::QCollection<1> oned_quadrature_collection; @@ -250,7 +286,8 @@ class DGBase virtual void assemble_volume_terms_implicit( const dealii::FEValues &fe_values_volume, const std::vector ¤t_dofs_indices, - dealii::Vector ¤t_cell_rhs) = 0; + dealii::Vector ¤t_cell_rhs, + const dealii::FEValues &fe_values_lagrange) = 0; /// Evaluate the integral over the cell edges that are on domain boundaries /** Compute both the right-hand side and the block of the Jacobian */ virtual void assemble_boundary_term_implicit( @@ -276,7 +313,8 @@ class DGBase virtual void assemble_volume_terms_explicit( const dealii::FEValues &fe_values_volume, const std::vector ¤t_dofs_indices, - dealii::Vector ¤t_cell_rhs) = 0; + dealii::Vector ¤t_cell_rhs, + const dealii::FEValues &fe_values_lagrange) = 0; /// Evaluate the integral over the cell edges that are on domain boundaries virtual void assemble_boundary_term_explicit( const unsigned int boundary_id, @@ -323,13 +361,6 @@ class DGBase // // const dealii::QGaussLobatto volume_quadrature; // // const dealii::QGaussLobatto face_quadrature; - /// Create mapping collection for initializer list - dealii::hp::MappingCollection create_mapping_collection(const unsigned int max_degree) const; - - /// Create FECollection for initializer list - dealii::hp::FECollection create_fe_collection(const unsigned int max_degree) const; - - /// Update flags needed at volume points. const dealii::UpdateFlags volume_update_flags = dealii::update_values | dealii::update_gradients | dealii::update_quadrature_points | dealii::update_JxW_values; /// Update flags needed at face points. @@ -338,26 +369,42 @@ class DGBase /** NOTE: With hp-adaptation, might need to query neighbor's quadrature points depending on the order of the cells. */ const dealii::UpdateFlags neighbor_face_update_flags = dealii::update_values | dealii::update_gradients; - dealii::hp::FEValues fe_values_collection_volume; ///< FEValues of volume. - dealii::hp::FEFaceValues fe_values_collection_face_int; ///< FEValues of interior face. - dealii::hp::FEFaceValues fe_values_collection_face_ext; ///< FEValues of exterior face. - dealii::hp::FESubfaceValues fe_values_collection_subface; ///< FEValues of subface. /// Lagrange basis used in strong form /** This is a collection of scalar Lagrange bases */ const dealii::hp::FECollection fe_collection_lagrange; - dealii::hp::FEValues fe_values_collection_volume_lagrange; + //dealii::hp::FEValues fe_values_collection_volume_lagrange; + + +public: + /// Degrees of freedom handler + /* Allows us to iterate over the finite elements' degrees of freedom. + * Note that since we are not using FESystem, we need to multiply + * the index by a factor of "nstate" + * + * Must be defined after fe_dg since it is a subscriptor of fe_dg. + * Destructor are called in reverse order in which they appear in class definition. + */ + dealii::hp::DoFHandler dof_handler; + + /// Degrees of freedom handler for the high-order grid + dealii::DoFHandler dof_handler_grid; + +protected: + /// Main underlying mapping responsible for all mapping evaluations + //dealii::MappingFEField, dealii::hp::DoFHandler > high_order_grid_mapping; + dealii::MappingFEField, dealii::DoFHandler > high_order_grid_mapping; MPI_Comm mpi_communicator; ///< MPI communicator dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 private: - /// Makes for cleaner doxygen documentation - using CollectionTuple = std::tuple< dealii::hp::MappingCollection, dealii::hp::FECollection, - dealii::hp::QCollection, dealii::hp::QCollection, dealii::hp::QCollection<1>, - dealii::hp::FECollection >; /// Used in the delegated constructor - CollectionTuple create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const; + /** The main reason we use this weird function is because all of the above objects + * need to be looped with the various p-orders. This function allows us to do this in a + * single function instead of having like 6 different functions to initialize each of them. + */ + MassiveCollectionTuple create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const; }; // end of DGBase class @@ -367,11 +414,17 @@ class DGBase template class DGWeak : public DGBase { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + using Triangulation = dealii::Triangulation; +#else + using Triangulation = dealii::parallel::distributed::Triangulation; +#endif public: /// Constructor DGWeak( const Parameters::AllParameters *const parameters_input, - const unsigned int degree); + const unsigned int degree, + Triangulation *const triangulation_input); ~DGWeak(); ///< Destructor @@ -397,7 +450,8 @@ class DGWeak : public DGBase void assemble_volume_terms_implicit( const dealii::FEValues &fe_values_volume, const std::vector ¤t_dofs_indices, - dealii::Vector ¤t_cell_rhs); + dealii::Vector ¤t_cell_rhs, + const dealii::FEValues &fe_values_lagrange); /// Evaluate the integral over the cell edges that are on domain boundaries /** Compute both the right-hand side and the block of the Jacobian */ void assemble_boundary_term_implicit( @@ -423,7 +477,8 @@ class DGWeak : public DGBase void assemble_volume_terms_explicit( const dealii::FEValues &fe_values_volume, const std::vector ¤t_dofs_indices, - dealii::Vector ¤t_cell_rhs); + dealii::Vector ¤t_cell_rhs, + const dealii::FEValues &fe_values_lagrange); /// Evaluate the integral over the cell edges that are on domain boundaries void assemble_boundary_term_explicit( const unsigned int boundary_id, @@ -452,11 +507,17 @@ class DGWeak : public DGBase template class DGStrong : public DGBase { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + using Triangulation = dealii::Triangulation; +#else + using Triangulation = dealii::parallel::distributed::Triangulation; +#endif public: /// Constructor DGStrong( const Parameters::AllParameters *const parameters_input, - const unsigned int degree); + const unsigned int degree, + Triangulation *const triangulation_input); /// Destructor ~DGStrong(); @@ -482,7 +543,8 @@ class DGStrong : public DGBase void assemble_volume_terms_implicit( const dealii::FEValues &fe_values_volume, const std::vector ¤t_dofs_indices, - dealii::Vector ¤t_cell_rhs); + dealii::Vector ¤t_cell_rhs, + const dealii::FEValues &fe_values_lagrange); /// Evaluate the integral over the cell edges that are on domain boundaries /** Compute both the right-hand side and the block of the Jacobian */ void assemble_boundary_term_implicit( @@ -508,7 +570,8 @@ class DGStrong : public DGBase void assemble_volume_terms_explicit( const dealii::FEValues &fe_values_volume, const std::vector ¤t_dofs_indices, - dealii::Vector ¤t_cell_rhs); + dealii::Vector ¤t_cell_rhs, + const dealii::FEValues &fe_values_lagrange); /// Evaluate the integral over the cell edges that are on domain boundaries void assemble_boundary_term_explicit( const unsigned int boundary_id, @@ -538,13 +601,19 @@ class DGStrong : public DGBase template class DGFactory { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + using Triangulation = dealii::Triangulation; +#else + using Triangulation = dealii::parallel::distributed::Triangulation; +#endif public: /// Creates a derived object DG, but returns it as DGBase. /** That way, the caller is agnostic to the number of state variables */ static std::shared_ptr< DGBase > create_discontinuous_galerkin( const Parameters::AllParameters *const parameters_input, - const unsigned int degree); + const unsigned int degree, + Triangulation *const triangulation_input); }; } // PHiLiP namespace diff --git a/src/dg/high_order_grid.cpp b/src/dg/high_order_grid.cpp new file mode 100644 index 000000000..2ede2e16d --- /dev/null +++ b/src/dg/high_order_grid.cpp @@ -0,0 +1,47 @@ +#include + +#include + +#include + +#include "high_order_grid.h" +namespace PHiLiP { + +template +HighOrderGrid::HighOrderGrid( + const Parameters::AllParameters *const parameters_input, + const unsigned int max_degree, + dealii::Triangulation *const triangulation_input) + : all_parameters(parameters_input) + , max_degree(max_degree) + , triangulation(triangulation_input) + , fe_q(max_degree) // The grid must be at least p1. A p0 solution required a p1 grid. + , fe_system(dealii::FESystem(fe_q,dim)) // The grid must be at least p1. A p0 solution required a p1 grid. + , mpi_communicator(MPI_COMM_WORLD) + , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) +{ + dof_handler_grid.initialize(*triangulation, fe_system); + dof_handler_grid.distribute_dofs(fe_system); + + locally_owned_dofs_grid = dof_handler_grid.locally_owned_dofs(); + dealii::DoFTools::extract_locally_relevant_dofs(dof_handler_grid, ghost_dofs_grid); + locally_relevant_dofs_grid = ghost_dofs_grid; + ghost_dofs_grid.subtract_set(locally_owned_dofs_grid); + nodes.reinit(locally_owned_dofs_grid, ghost_dofs_grid, mpi_communicator); + +} + +template +dealii::MappingFEField +HighOrderGrid::get_MappingFEField() { + const dealii::ComponentMask mask(dim, true); + dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask); + + dealii::MappingFEField mapping(dof_handler_grid,nodes,mask); + + return mapping; +} + +//template class HighOrderGrid; +template class HighOrderGrid, dealii::DoFHandler>; +} // namespace PHiLiP diff --git a/src/dg/high_order_grid.h b/src/dg/high_order_grid.h new file mode 100644 index 000000000..1b478a7d2 --- /dev/null +++ b/src/dg/high_order_grid.h @@ -0,0 +1,71 @@ +#ifndef __HIGHORDERGRID_H__ +#define __HIGHORDERGRID_H__ + +#include + +#include +#include + +#include + +#include + +#include "parameters/all_parameters.h" + +namespace PHiLiP { + +/** This class is used to generate and control the high-order nodes given a Triangulation and a Manifold. + * This will especially be useful when performing shape optimization, where the surface and volume + * nodes will need to be displaced. + */ +//template +template , typename DoFHandlerType = dealii::DoFHandler> +class HighOrderGrid +{ +public: + /// Principal constructor that will call delegated constructor. + HighOrderGrid(const Parameters::AllParameters *const parameters_input, const unsigned int max_degree, dealii::Triangulation *const triangulation_input); + + dealii::MappingFEField get_MappingFEField(); + + const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters + + /// Maximum degree of the grid. + const unsigned int max_degree; + + dealii::Triangulation *triangulation; ///< Mesh + + /// Get evaluate high-order grid from Triangulation and Manifold + + + /// Degrees of freedom handler for the high-order grid + dealii::DoFHandler dof_handler_grid; + + /// Current nodal coefficients of the high-order grid. + dealii::LinearAlgebra::distributed::Vector nodes; + + /// List of surface points + dealii::LinearAlgebra::distributed::Vector surface_nodes; + + + /// RBF mesh deformation + void deform_mesh(dealii::LinearAlgebra::distributed::Vector surface_displacements); + + + /// Using system of FE_Q to represent the grid + const dealii::FE_Q fe_q; + const dealii::FESystem fe_system; + + dealii::IndexSet locally_owned_dofs_grid; ///< Locally own degrees of freedom for the grid + dealii::IndexSet ghost_dofs_grid; ///< Locally relevant ghost degrees of freedom for the grid + dealii::IndexSet locally_relevant_dofs_grid; ///< Union of locally owned degrees of freedom and relevant ghost degrees of freedom for the grid +protected: + + MPI_Comm mpi_communicator; ///< MPI communicator + dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 + +}; + +} // namespace PHiLiP + +#endif diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp index 785a8cad8..a3ae702fc 100644 --- a/src/dg/strong_dg.cpp +++ b/src/dg/strong_dg.cpp @@ -19,12 +19,19 @@ namespace PHiLiP { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + template using Triangulation = dealii::Triangulation; +#else + template using Triangulation = dealii::parallel::distributed::Triangulation; +#endif + template DGStrong::DGStrong( const Parameters::AllParameters *const parameters_input, - const unsigned int degree) - : DGBase::DGBase(nstate, parameters_input, degree) // Use DGBase constructor + const unsigned int degree, + Triangulation *const triangulation_input) + : DGBase::DGBase(nstate, parameters_input, degree, triangulation_input) // Use DGBase constructor { using ADtype = Sacado::Fad::DFad; pde_physics = Physics::PhysicsFactory ::create_Physics(parameters_input); @@ -48,7 +55,8 @@ template void DGStrong::assemble_volume_terms_implicit( const dealii::FEValues &fe_values_vol, const std::vector &cell_dofs_indices, - dealii::Vector &local_rhs_int_cell) + dealii::Vector &local_rhs_int_cell, + const dealii::FEValues &fe_values_lagrange) { using ADtype = Sacado::Fad::DFad; @@ -108,7 +116,7 @@ void DGStrong::assemble_volume_terms_implicit( // Evaluate flux divergence by interpolating the flux // Since we have nodal values of the flux, we use the Lagrange polynomials to obtain the gradients at the quadrature points. - const dealii::FEValues &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values(); + //const dealii::FEValues &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values(); std::vector flux_divergence(n_quad_pts); std::array,nstate>,dim> f; @@ -477,7 +485,8 @@ template void DGStrong::assemble_volume_terms_explicit( const dealii::FEValues &fe_values_vol, const std::vector &cell_dofs_indices, - dealii::Vector &local_rhs_int_cell) + dealii::Vector &local_rhs_int_cell, + const dealii::FEValues &fe_values_lagrange) { //std::cout << "assembling cell terms" << std::endl; using ADtype = Sacado::Fad::DFad; @@ -536,7 +545,7 @@ void DGStrong::assemble_volume_terms_explicit( // Evaluate flux divergence by interpolating the flux // Since we have nodal values of the flux, we use the Lagrange polynomials to obtain the gradients at the quadrature points. - const dealii::FEValues &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values(); + //const dealii::FEValues &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values(); std::vector flux_divergence(n_quad_pts); for (int istate = 0; istate does not work for 1D + template using Triangulation = dealii::Triangulation; +#else + template using Triangulation = dealii::parallel::distributed::Triangulation; +#endif + template DGWeak::DGWeak( const Parameters::AllParameters *const parameters_input, - const unsigned int degree) - : DGBase::DGBase(nstate, parameters_input, degree) // Use DGBase constructor + const unsigned int degree, + Triangulation *const triangulation_input) + : DGBase::DGBase(nstate, parameters_input, degree, triangulation_input) // Use DGBase constructor { using ADtype = Sacado::Fad::DFad; pde_physics = Physics::PhysicsFactory ::create_Physics(parameters_input); @@ -45,7 +52,8 @@ template void DGWeak::assemble_volume_terms_implicit( const dealii::FEValues &fe_values_vol, const std::vector &cell_dofs_indices, - dealii::Vector &local_rhs_int_cell) + dealii::Vector &local_rhs_int_cell, + const dealii::FEValues &/*fe_values_lagrange*/) { using ADtype = Sacado::Fad::DFad; using ADArray = std::array; @@ -431,7 +439,8 @@ template void DGWeak::assemble_volume_terms_explicit( const dealii::FEValues &fe_values_vol, const std::vector &cell_dofs_indices, - dealii::Vector &local_rhs_int_cell) + dealii::Vector &local_rhs_int_cell, + const dealii::FEValues &/*fe_values_lagrange*/) { using doubleArray = std::array; using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,real>, nstate >; diff --git a/src/testing/burgers_stability.cpp b/src/testing/burgers_stability.cpp index 5b2e459a3..56d1516e0 100644 --- a/src/testing/burgers_stability.cpp +++ b/src/testing/burgers_stability.cpp @@ -55,9 +55,8 @@ int BurgersEnergyStability::run_test() const dealii::GridGenerator::hyper_cube(grid, left, right, colorize); grid.refine_global(n_refinements); pcout << "Grid generated and refined" << std::endl; - std::shared_ptr < PHiLiP::DGBase > dg = PHiLiP::DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree); + std::shared_ptr < PHiLiP::DGBase > dg = PHiLiP::DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree, &grid); pcout << "dg created" <set_triangulation(&grid); dg->allocate_system (); pcout << "Implement initial conditions" << std::endl; diff --git a/src/testing/euler_cylinder.cpp b/src/testing/euler_cylinder.cpp index 7a003294f..be96398ec 100644 --- a/src/testing/euler_cylinder.cpp +++ b/src/testing/euler_cylinder.cpp @@ -189,8 +189,7 @@ ::run_test () const } // Create DG object - std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree, &grid); // Initialize coarse grid solution with free-stream diff --git a/src/testing/euler_entropy_waves.cpp b/src/testing/euler_entropy_waves.cpp index 31ac2397c..ec84483c9 100644 --- a/src/testing/euler_entropy_waves.cpp +++ b/src/testing/euler_entropy_waves.cpp @@ -188,8 +188,7 @@ ::run_test () const if (random_factor > 0.0) dealii::GridTools::distort_random (random_factor, grid, keep_boundary); // Create DG object using the factory - std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree, &grid); dg->allocate_system (); // Initialize solution with vortex function at time t=0 diff --git a/src/testing/euler_gaussian_bump.cpp b/src/testing/euler_gaussian_bump.cpp index cc1a0371b..1bb4ecb20 100644 --- a/src/testing/euler_gaussian_bump.cpp +++ b/src/testing/euler_gaussian_bump.cpp @@ -222,8 +222,7 @@ ::run_test () const grid.set_manifold ( manifold_id, bump_manifold ); // Create DG object - std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree, &grid); // Initialize coarse grid solution with free-stream diff --git a/src/testing/euler_vortex.cpp b/src/testing/euler_vortex.cpp index 4a3b5ccba..29624e2e2 100644 --- a/src/testing/euler_vortex.cpp +++ b/src/testing/euler_vortex.cpp @@ -230,8 +230,7 @@ ::run_test () const if (random_factor > 0.0) dealii::GridTools::distort_random (random_factor, grid, keep_boundary); // Create DG object using the factory - std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree, &grid); dg->allocate_system (); // Initialize solution with vortex function at time t=0 diff --git a/src/testing/grid_study.cpp b/src/testing/grid_study.cpp index 84d6f7539..ecd46844c 100644 --- a/src/testing/grid_study.cpp +++ b/src/testing/grid_study.cpp @@ -141,8 +141,7 @@ ::run_test () const dealii::Triangulation::smoothing_on_coarsening)); #endif dealii::GridGenerator::subdivided_hyper_cube(grid_super_fine, n_1d_cells[n_grids_input-1]); - std::shared_ptr < DGBase > dg_super_fine = DGFactory::create_discontinuous_galerkin(¶m, p_end); - dg_super_fine->set_triangulation(&grid_super_fine); + std::shared_ptr < DGBase > dg_super_fine = DGFactory::create_discontinuous_galerkin(¶m, p_end, &grid_super_fine); dg_super_fine->allocate_system (); initialize_perturbed_solution(*dg_super_fine, *physics_double); @@ -233,8 +232,7 @@ ::run_test () const using ADtype = Sacado::Fad::DFad; // Create DG object using the factory - std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree, &grid); dg->allocate_system (); //dg->evaluate_inverse_mass_matrices(); // diff --git a/tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp b/tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp index 69832e0b4..c307f8232 100644 --- a/tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp +++ b/tests/integration_tests_control_files/advection_explicit_periodic/advection_explicit_periodic.cpp @@ -72,8 +72,7 @@ int AdvectionPeriodic::run_test() grid.refine_global(n_refinements); - std::shared_ptr < PHiLiP::DGBase > dg = PHiLiP::DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < PHiLiP::DGBase > dg = PHiLiP::DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree, &grid); dg->allocate_system (); // for (auto current_cell = dg->dof_handler.begin_active(); current_cell != dg->dof_handler.end(); ++current_cell) { diff --git a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp b/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp index 9d104faee..0f3fd1d64 100644 --- a/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp +++ b/tests/integration_tests_control_files/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp @@ -143,8 +143,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); - dg->set_triangulation(&grid); + std::shared_ptr < PHiLiP::DGBase > dg = PHiLiP::DGFactory::create_discontinuous_galerkin(all_parameters, poly_degree, &grid); dg->allocate_system (); std::cout << "Implement initial conditions" << std::endl; diff --git a/tests/unit_tests/CMakeLists.txt b/tests/unit_tests/CMakeLists.txt index bd4057f01..8b7780df4 100644 --- a/tests/unit_tests/CMakeLists.txt +++ b/tests/unit_tests/CMakeLists.txt @@ -2,3 +2,4 @@ add_subdirectory(numerical_flux) add_subdirectory(regression) add_subdirectory(euler_unit_test) +add_subdirectory(grid) diff --git a/tests/unit_tests/grid/CMakeLists.txt b/tests/unit_tests/grid/CMakeLists.txt new file mode 100644 index 000000000..287d37dfb --- /dev/null +++ b/tests/unit_tests/grid/CMakeLists.txt @@ -0,0 +1,36 @@ +set(TEST_SRC + high_order_grid_test.cpp + ) + +foreach(dim RANGE 2 3) + + # Output executable + string(CONCAT TEST_TARGET ${dim}D_HighOrder_MappingFEField) + 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(ParametersLib ParametersLibrary) + string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) + target_link_libraries(${TEST_TARGET} ${ParametersLib}) + target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) + # Setup target with deal.II + DEAL_II_SETUP_TARGET(${TEST_TARGET}) + + add_test( + NAME ${TEST_TARGET} + COMMAND mpirun -n 1 ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} + ) + + unset(TEST_TARGET) + unset(DiscontinuousGalerkinLib) + unset(ParametersLib) + +endforeach() diff --git a/tests/unit_tests/grid/high_order_grid_test.cpp b/tests/unit_tests/grid/high_order_grid_test.cpp new file mode 100644 index 000000000..313249216 --- /dev/null +++ b/tests/unit_tests/grid/high_order_grid_test.cpp @@ -0,0 +1,227 @@ + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include + + +#include +#include + +#include "dg/high_order_grid.h" +#include "parameters/all_parameters.h" + +// https://en.wikipedia.org/wiki/Volume_of_an_n-ball#Recursions +template +double volume_n_ball(const double radius) +{ + const double pi = dealii::numbers::PI; + return (2.0*pi*radius*radius/dim)*volume_n_ball(radius); +} +template <> double volume_n_ball<0>(const double /*radius*/) { return 1.0; } +template <> double volume_n_ball<1>(const double radius) { return 2.0*radius; } + +int main (int argc, char * argv[]) +{ + const int dim = PHILIP_DIM; + int success_bool = true; + + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); + + using namespace PHiLiP; + + dealii::ParameterHandler parameter_handler; + Parameters::AllParameters::declare_parameters (parameter_handler); + Parameters::AllParameters all_parameters; + all_parameters.parse_parameters (parameter_handler); + + std::vector fail_conv_poly; + std::vector fail_conv_slop; + std::vector convergence_table_vector; + + const unsigned int p_start = 1; + const unsigned int p_end = 4; + const unsigned int n_grids = 3; + //const std::vector n_1d_cells = {2,4,8,16}; + + //const unsigned int n_cells_circle = n_1d_cells[0]; + //const unsigned int n_cells_radial = 3*n_cells_circle; + + const dealii::Point center; // Constructor initializes with 0; + const double inner_radius = 1, outer_radius = inner_radius*10; + const double exact_volume = (volume_n_ball(outer_radius) - volume_n_ball(inner_radius)) / std::pow(2,dim); + + for (unsigned int poly_degree = p_start; poly_degree <= p_end; ++poly_degree) { + + std::vector grid_size(n_grids); + + dealii::ConvergenceTable convergence_table; + + // Generate grid and mapping + dealii::parallel::distributed::Triangulation grid( + MPI_COMM_WORLD, + typename dealii::Triangulation::MeshSmoothing( + dealii::Triangulation::smoothing_on_refinement | + dealii::Triangulation::smoothing_on_coarsening)); + + //dealii::Triangulation grid; + + const int n_cells = 0*(dim-1); + dealii::GridGenerator::quarter_hyper_shell(grid, center, inner_radius, outer_radius, n_cells);//, n_cells = 0, colorize = false); + //dealii::GridGenerator::hyper_shell(grid, center, inner_radius, outer_radius, n_cells);//, n_cells = 0, colorize = false); + grid.set_all_manifold_ids(0); + grid.set_manifold(0, dealii::SphericalManifold(center)); + + std::vector volume_error(n_grids); + + for (unsigned int igrid=0; igrid high_order_grid(&all_parameters, poly_degree, &grid); + + const unsigned int n_dofs = high_order_grid.dof_handler_grid.n_dofs(); + + const unsigned int n_global_active_cells = grid.n_global_active_cells(); + + // Overintegrate the error to make sure there is not integration error in the error estimate + const int overintegrate = 0; + const unsigned int n_quad_pts_1D = poly_degree+1+overintegrate; + dealii::QGauss quadrature(n_quad_pts_1D); + const unsigned int n_quad_pts = quadrature.size(); + //const bool use_mapping_q_on_all_cells = true; + //dealii::MappingQ mapping(poly_degree, use_mapping_q_on_all_cells); + //dealii::Mapping mapping = high_order_grid.get_MappingFEField(); + dealii::MappingFEField, dealii::DoFHandler> mapping = high_order_grid.get_MappingFEField(); + + //grid.reset_all_manifolds(); + + // Integrate solution error and output error + dealii::FEValues fe_values(mapping, high_order_grid.fe_system, quadrature, + dealii::update_jacobians + //| dealii::update_volume_elements + | dealii::update_JxW_values + | dealii::update_quadrature_points); + + std::ofstream out_before("before_move_grid-" + std::to_string(poly_degree) + "-" + std::to_string(igrid) + ".vtk"); + dealii::GridOut grid_out_before; + grid_out_before.write_vtk(grid, out_before); + for (auto cell = high_order_grid.dof_handler_grid.begin_active(); cell!=high_order_grid.dof_handler_grid.end(); ++cell) { + + if (!cell->is_locally_owned()) continue; + const unsigned int fe_index_cell = cell->active_fe_index(); + const dealii::FESystem &fe_ref = high_order_grid.fe_system[fe_index_cell]; + const unsigned int n_dofs_cell = fe_ref.n_dofs_per_cell(); + + // Obtain the mapping from local dof indices to global dof indices + std::vector dofs_indices; + dofs_indices.resize(n_dofs_cell); + cell->get_dof_indices (dofs_indices); + + for (unsigned int idof=0; idofis_locally_owned()) continue; + n_cell++; + + fe_values.reinit (cell); + + double cell_volume = 0.0; + for (unsigned int iquad=0; iquad 2 ) { + before_last_slope = log(volume_error[n_grids-2]/volume_error[n_grids-3]) + / log(grid_size[n_grids-2]/grid_size[n_grids-3]); + } + const double slope_avg = 0.5*(before_last_slope+last_slope); + const double slope_diff = slope_avg-expected_slope; + + const double slope_deficit_tolerance = -0.1;//-std::abs(manu_grid_conv_param.slope_deficit_tolerance); + + if (slope_diff < slope_deficit_tolerance) { + pcout << std::endl << "Convergence order not achieved. Average last 2 slopes of " << slope_avg << " instead of expected " + << expected_slope << " within a tolerance of " << slope_deficit_tolerance << std::endl; + if(poly_degree!=0) fail_conv_poly.push_back(poly_degree); + if(poly_degree!=0) fail_conv_slop.push_back(slope_avg); + } + + } + pcout << std::endl << std::endl << std::endl << std::endl; + pcout << " ********************************************" << std::endl; + pcout << " Convergence summary" << std::endl; + pcout << " ********************************************" << std::endl; + for (auto conv = convergence_table_vector.begin(); conv!=convergence_table_vector.end(); conv++) { + if (pcout.is_active()) conv->write_text(pcout.get_stream()); + pcout << " ********************************************" << std::endl; + } + int n_fail_poly = fail_conv_poly.size(); + if (n_fail_poly > 0) { + for (int ifail=0; ifail < n_fail_poly; ++ifail) { + const double expected_slope = fail_conv_poly[ifail]+1; + const double slope_deficit_tolerance = -0.1; + pcout << std::endl << "Convergence order not achieved for polynomial p = " << fail_conv_poly[ifail] + << ". Slope of " << fail_conv_slop[ifail] << " instead of expected " << expected_slope + << " within a tolerance of " << slope_deficit_tolerance << std::endl; + } + + success_bool = false; + } + return success_bool; +} diff --git a/tests/unit_tests/regression/jacobian_matrix_regression.cpp b/tests/unit_tests/regression/jacobian_matrix_regression.cpp index 50dc1a26c..5e357f53b 100644 --- a/tests/unit_tests/regression/jacobian_matrix_regression.cpp +++ b/tests/unit_tests/regression/jacobian_matrix_regression.cpp @@ -62,8 +62,7 @@ int main (int argc, char * argv[]) // Assemble Jacobian all_parameters.pde_type = *pde; - std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(&all_parameters, poly_degree); - dg->set_triangulation(&grid); + std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(&all_parameters, poly_degree, &grid); dg->allocate_system (); dg->assemble_residual(true);