Skip to content

Commit

Permalink
Add basic HighOrderGrid class and test.
Browse files Browse the repository at this point in the history
  • Loading branch information
dougshidong committed Sep 23, 2019
1 parent 38ce54d commit 4dd6d62
Show file tree
Hide file tree
Showing 20 changed files with 812 additions and 145 deletions.
14 changes: 11 additions & 3 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/dg/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
set(DG_SOURCE
high_order_grid.cpp
dg.cpp
weak_dg.cpp
strong_dg.cpp
Expand Down
317 changes: 258 additions & 59 deletions src/dg/dg.cpp

Large diffs are not rendered by default.

175 changes: 122 additions & 53 deletions src/dg/dg.h

Large diffs are not rendered by default.

47 changes: 47 additions & 0 deletions src/dg/high_order_grid.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
#include <deal.II/fe/fe_q.h>

#include <deal.II/dofs/dof_tools.h>

#include <deal.II/numerics/vector_tools.h>

#include "high_order_grid.h"
namespace PHiLiP {

template <int dim, typename real, typename VectorType , typename DoFHandlerType>
HighOrderGrid<dim,real,VectorType,DoFHandlerType>::HighOrderGrid(
const Parameters::AllParameters *const parameters_input,
const unsigned int max_degree,
dealii::Triangulation<dim> *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<dim>(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 <int dim, typename real, typename VectorType , typename DoFHandlerType>
dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>
HighOrderGrid<dim,real,VectorType,DoFHandlerType>::get_MappingFEField() {
const dealii::ComponentMask mask(dim, true);
dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask);

dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> mapping(dof_handler_grid,nodes,mask);

return mapping;
}

//template class HighOrderGrid<PHILIP_DIM, double>;
template class HighOrderGrid<PHILIP_DIM, double, dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<PHILIP_DIM>>;
} // namespace PHiLiP
71 changes: 71 additions & 0 deletions src/dg/high_order_grid.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
#ifndef __HIGHORDERGRID_H__
#define __HIGHORDERGRID_H__

#include <deal.II/grid/tria.h>

#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_fe_field.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/lac/vector.h>

#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 <int dim, typename real, typename VectorType , typename DoFHandlerType>
template <int dim = PHILIP_DIM, typename real = double, typename VectorType = dealii::LinearAlgebra::distributed::Vector<double>, typename DoFHandlerType = dealii::DoFHandler<PHILIP_DIM>>
class HighOrderGrid
{
public:
/// Principal constructor that will call delegated constructor.
HighOrderGrid(const Parameters::AllParameters *const parameters_input, const unsigned int max_degree, dealii::Triangulation<dim> *const triangulation_input);

dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> get_MappingFEField();

const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters

/// Maximum degree of the grid.
const unsigned int max_degree;

dealii::Triangulation<dim> *triangulation; ///< Mesh

/// Get evaluate high-order grid from Triangulation and Manifold


/// Degrees of freedom handler for the high-order grid
dealii::DoFHandler<dim> dof_handler_grid;

/// Current nodal coefficients of the high-order grid.
dealii::LinearAlgebra::distributed::Vector<double> nodes;

/// List of surface points
dealii::LinearAlgebra::distributed::Vector<double> surface_nodes;


/// RBF mesh deformation
void deform_mesh(dealii::LinearAlgebra::distributed::Vector<double> surface_displacements);


/// Using system of FE_Q to represent the grid
const dealii::FE_Q<dim> fe_q;
const dealii::FESystem<dim> 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
21 changes: 15 additions & 6 deletions src/dg/strong_dg.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,19 @@

namespace PHiLiP {

#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
template <int dim> using Triangulation = dealii::Triangulation<dim>;
#else
template <int dim> using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
#endif


template <int dim, int nstate, typename real>
DGStrong<dim,nstate,real>::DGStrong(
const Parameters::AllParameters *const parameters_input,
const unsigned int degree)
: DGBase<dim,real>::DGBase(nstate, parameters_input, degree) // Use DGBase constructor
const unsigned int degree,
Triangulation *const triangulation_input)
: DGBase<dim,real>::DGBase(nstate, parameters_input, degree, triangulation_input) // Use DGBase constructor
{
using ADtype = Sacado::Fad::DFad<real>;
pde_physics = Physics::PhysicsFactory<dim,nstate,ADtype> ::create_Physics(parameters_input);
Expand All @@ -48,7 +55,8 @@ template <int dim, int nstate, typename real>
void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
const dealii::FEValues<dim,dim> &fe_values_vol,
const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
dealii::Vector<real> &local_rhs_int_cell)
dealii::Vector<real> &local_rhs_int_cell,
const dealii::FEValues<dim,dim> &fe_values_lagrange)
{

using ADtype = Sacado::Fad::DFad<real>;
Expand Down Expand Up @@ -108,7 +116,7 @@ void DGStrong<dim,nstate,real>::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<dim,dim> &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values();
//const dealii::FEValues<dim,dim> &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values();
std::vector<ADArray> flux_divergence(n_quad_pts);

std::array<std::array<std::vector<ADtype>,nstate>,dim> f;
Expand Down Expand Up @@ -477,7 +485,8 @@ template <int dim, int nstate, typename real>
void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
const dealii::FEValues<dim,dim> &fe_values_vol,
const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
dealii::Vector<real> &local_rhs_int_cell)
dealii::Vector<real> &local_rhs_int_cell,
const dealii::FEValues<dim,dim> &fe_values_lagrange)
{
//std::cout << "assembling cell terms" << std::endl;
using ADtype = Sacado::Fad::DFad<real>;
Expand Down Expand Up @@ -536,7 +545,7 @@ void DGStrong<dim,nstate,real>::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<dim,dim> &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values();
//const dealii::FEValues<dim,dim> &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values();
std::vector<ADArray> flux_divergence(n_quad_pts);
for (int istate = 0; istate<nstate; ++istate) {
for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
Expand Down
17 changes: 13 additions & 4 deletions src/dg/weak_dg.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,11 +17,18 @@

namespace PHiLiP {

#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
template <int dim> using Triangulation = dealii::Triangulation<dim>;
#else
template <int dim> using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
#endif

template <int dim, int nstate, typename real>
DGWeak<dim,nstate,real>::DGWeak(
const Parameters::AllParameters *const parameters_input,
const unsigned int degree)
: DGBase<dim,real>::DGBase(nstate, parameters_input, degree) // Use DGBase constructor
const unsigned int degree,
Triangulation *const triangulation_input)
: DGBase<dim,real>::DGBase(nstate, parameters_input, degree, triangulation_input) // Use DGBase constructor
{
using ADtype = Sacado::Fad::DFad<real>;
pde_physics = Physics::PhysicsFactory<dim,nstate,ADtype> ::create_Physics(parameters_input);
Expand All @@ -45,7 +52,8 @@ template <int dim, int nstate, typename real>
void DGWeak<dim,nstate,real>::assemble_volume_terms_implicit(
const dealii::FEValues<dim,dim> &fe_values_vol,
const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
dealii::Vector<real> &local_rhs_int_cell)
dealii::Vector<real> &local_rhs_int_cell,
const dealii::FEValues<dim,dim> &/*fe_values_lagrange*/)
{
using ADtype = Sacado::Fad::DFad<real>;
using ADArray = std::array<ADtype,nstate>;
Expand Down Expand Up @@ -431,7 +439,8 @@ template <int dim, int nstate, typename real>
void DGWeak<dim,nstate,real>::assemble_volume_terms_explicit(
const dealii::FEValues<dim,dim> &fe_values_vol,
const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
dealii::Vector<real> &local_rhs_int_cell)
dealii::Vector<real> &local_rhs_int_cell,
const dealii::FEValues<dim,dim> &/*fe_values_lagrange*/)
{
using doubleArray = std::array<real,nstate>;
using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,real>, nstate >;
Expand Down
3 changes: 1 addition & 2 deletions src/testing/burgers_stability.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,8 @@ int BurgersEnergyStability<dim, nstate>::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<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree);
std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree, &grid);
pcout << "dg created" <<std::endl;
dg->set_triangulation(&grid);
dg->allocate_system ();

pcout << "Implement initial conditions" << std::endl;
Expand Down
3 changes: 1 addition & 2 deletions src/testing/euler_cylinder.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -189,8 +189,7 @@ ::run_test () const
}

// Create DG object
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, &grid);


// Initialize coarse grid solution with free-stream
Expand Down
3 changes: 1 addition & 2 deletions src/testing/euler_entropy_waves.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, &grid);
dg->allocate_system ();

// Initialize solution with vortex function at time t=0
Expand Down
3 changes: 1 addition & 2 deletions src/testing/euler_gaussian_bump.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -222,8 +222,7 @@ ::run_test () const
grid.set_manifold ( manifold_id, bump_manifold );

// Create DG object
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, &grid);


// Initialize coarse grid solution with free-stream
Expand Down
3 changes: 1 addition & 2 deletions src/testing/euler_vortex.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, &grid);
dg->allocate_system ();

// Initialize solution with vortex function at time t=0
Expand Down
6 changes: 2 additions & 4 deletions src/testing/grid_study.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -141,8 +141,7 @@ ::run_test () const
dealii::Triangulation<dim>::smoothing_on_coarsening));
#endif
dealii::GridGenerator::subdivided_hyper_cube(grid_super_fine, n_1d_cells[n_grids_input-1]);
std::shared_ptr < DGBase<dim, double> > dg_super_fine = DGFactory<dim,double>::create_discontinuous_galerkin(&param, p_end);
dg_super_fine->set_triangulation(&grid_super_fine);
std::shared_ptr < DGBase<dim, double> > dg_super_fine = DGFactory<dim,double>::create_discontinuous_galerkin(&param, p_end, &grid_super_fine);
dg_super_fine->allocate_system ();

initialize_perturbed_solution(*dg_super_fine, *physics_double);
Expand Down Expand Up @@ -233,8 +232,7 @@ ::run_test () const
using ADtype = Sacado::Fad::DFad<double>;

// Create DG object using the factory
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, &grid);
dg->allocate_system ();
//dg->evaluate_inverse_mass_matrices();
//
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -72,8 +72,7 @@ int AdvectionPeriodic<dim, nstate>::run_test()

grid.refine_global(n_refinements);

std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::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) {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -143,8 +143,7 @@ int EulerTaylorGreen<dim, nstate>::run_test()

grid.refine_global(n_refinements);

std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree);
dg->set_triangulation(&grid);
std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree, &grid);
dg->allocate_system ();

std::cout << "Implement initial conditions" << std::endl;
Expand Down
1 change: 1 addition & 0 deletions tests/unit_tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,4 @@
add_subdirectory(numerical_flux)
add_subdirectory(regression)
add_subdirectory(euler_unit_test)
add_subdirectory(grid)
36 changes: 36 additions & 0 deletions tests/unit_tests/grid/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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()
Loading

0 comments on commit 4dd6d62

Please sign in to comment.