[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
surface_GCL_test.cpp
1 #include <iomanip>
2 #include <cmath>
3 #include <limits>
4 #include <type_traits>
5 #include <math.h>
6 #include <iostream>
7 #include <stdlib.h>
8 
9 #include <deal.II/distributed/solution_transfer.h>
10 
11 #include "testing/tests.h"
12 
13 #include<fstream>
14 #include <deal.II/base/parameter_handler.h>
15 #include <deal.II/base/tensor.h>
16 #include <deal.II/numerics/vector_tools.h>
17 
18 #include <deal.II/grid/grid_generator.h>
19 #include <deal.II/grid/grid_refinement.h>
20 #include <deal.II/grid/grid_tools.h>
21 #include <deal.II/grid/grid_out.h>
22 #include <deal.II/grid/grid_in.h>
23 
24 #include <deal.II/dofs/dof_handler.h>
25 #include <deal.II/dofs/dof_tools.h>
26 #include <deal.II/dofs/dof_renumbering.h>
27 
28 #include <deal.II/dofs/dof_accessor.h>
29 
30 #include <deal.II/lac/vector.h>
31 #include <deal.II/lac/dynamic_sparsity_pattern.h>
32 #include <deal.II/lac/sparse_matrix.h>
33 
34 #include <deal.II/meshworker/dof_info.h>
35 
36 // Finally, we take our exact solution from the library as well as volume_quadrature
37 // and additional tools.
38 #include <deal.II/fe/mapping_q.h>
39 #include <deal.II/grid/manifold_lib.h>
40 #include <deal.II/numerics/data_out.h>
41 #include <deal.II/numerics/data_out_dof_data.h>
42 #include <deal.II/numerics/vector_tools.h>
43 #include <deal.II/numerics/vector_tools.templates.h>
44 
45 #include "dg/dg_base.hpp"
46 #include "dg/dg_factory.hpp"
47 #include "operators/operators.h"
48 #include "parameters/all_parameters.h"
49 #include "parameters/parameters.h"
50 //#include <GCL_test.h>
51 
52 const double TOLERANCE = 1E-6;
53 using namespace std;
54 //namespace PHiLiP {
55 
56 template <int dim>
57 class CurvManifold: public dealii::ChartManifold<dim,dim,dim> {
58  virtual dealii::Point<dim> pull_back(const dealii::Point<dim> &space_point) const override;
59  virtual dealii::Point<dim> push_forward(const dealii::Point<dim> &chart_point) const override;
60  virtual dealii::DerivativeForm<1,dim,dim> push_forward_gradient(const dealii::Point<dim> &chart_point) const override;
61 
62  virtual std::unique_ptr<dealii::Manifold<dim,dim> > clone() const override;
63 };
64 
65 template<int dim>
66 dealii::Point<dim> CurvManifold<dim>::pull_back(const dealii::Point<dim> &space_point) const
67 {
68  using namespace PHiLiP;
69  const double pi = atan(1)*4.0;
70  dealii::Point<dim> x_ref;
71  dealii::Point<dim> x_phys;
72  for(int idim=0; idim<dim; idim++){
73  x_ref[idim] = space_point[idim];
74  x_phys[idim] = space_point[idim];
75  }
76  dealii::Vector<double> function(dim);
77  dealii::FullMatrix<double> derivative(dim);
78  double beta =1.0/10.0;
79  double alpha =1.0/10.0;
80  int flag =0;
81  while(flag != dim){
82  if(dim==2){
83  function[0] = x_ref[0] - x_phys[0] +beta*std::cos(pi/2.0*x_ref[0])*std::cos(3.0*pi/2.0*x_ref[1]);
84  function[1] = x_ref[1] - x_phys[1] +beta*std::sin(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
85  }
86  else{
87  function[0] = x_ref[0] - x_phys[0] +alpha*(std::cos(pi * x_ref[2]) + std::cos(pi * x_ref[1]));
88  function[1] = x_ref[1] - x_phys[1] +alpha*exp(1.0-x_ref[1])*(std::sin(pi * x_ref[0]) + std::sin(pi* x_ref[2]));
89  function[2] = x_ref[2] - x_phys[2] +1.0/20.0*( std::sin(2.0 * pi * x_ref[0]) + std::sin(2.0 * pi * x_ref[1]));
90  }
91 
92  if(dim==2){
93  derivative[0][0] = 1.0 - beta* pi/2.0 * std::sin(pi/2.0*x_ref[0])*std::cos(3.0*pi/2.0*x_ref[1]);
94  derivative[0][1] = - beta*3.0 *pi/2.0 * std::cos(pi/2.0*x_ref[0])*std::sin(3.0*pi/2.0*x_ref[1]);
95 
96  derivative[1][0] = beta*2.0*pi*std::cos(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
97  derivative[1][1] = 1.0 -beta*pi/2.0*std::sin(2.0*pi*(x_ref[0]))*std::sin(pi/2.0*x_ref[1]);
98  }
99  else{
100  derivative[0][0] = 1.0;
101  derivative[0][1] = - alpha*pi*std::sin(pi*x_ref[1]);
102  derivative[0][2] = - alpha*pi*std::sin(pi*x_ref[2]);
103 
104  derivative[1][0] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[0]);
105  derivative[1][1] = 1.0 -alpha*exp(1.0-x_ref[1])*(std::sin(pi*x_ref[0])+std::sin(pi*x_ref[2]));
106  derivative[1][2] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[2]);
107  derivative[2][0] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[0]);
108  derivative[2][1] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[1]);
109  derivative[2][2] = 1.0;
110  }
111 
112  dealii::FullMatrix<double> Jacobian_inv(dim);
113  Jacobian_inv.invert(derivative);
114  dealii::Vector<double> Newton_Step(dim);
115  Jacobian_inv.vmult(Newton_Step, function);
116  for(int idim=0; idim<dim; idim++){
117  x_ref[idim] -= Newton_Step[idim];
118  }
119  flag=0;
120  for(int idim=0; idim<dim; idim++){
121  if(std::abs(function[idim]) < 1e-15)
122  flag++;
123  }
124  if(flag == dim)
125  break;
126  }
127  std::vector<double> function_check(dim);
128  if(dim==2){
129  function_check[0] = x_ref[0] + beta*std::cos(pi/2.0*x_ref[0])*std::cos(3.0*pi/2.0*x_ref[1]);
130  function_check[1] = x_ref[1] + beta*std::sin(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
131  }
132  else{
133  function_check[0] = x_ref[0] +alpha*(std::cos(pi * x_ref[2]) + std::cos(pi * x_ref[1]));
134  function_check[1] = x_ref[1] +alpha*exp(1.0-x_ref[1])*(std::sin(pi * x_ref[0]) + std::sin(pi* x_ref[2]));
135  function_check[2] = x_ref[2] +1.0/20.0*( std::sin(2.0 * pi * x_ref[0]) + std::sin(2.0 * pi * x_ref[1]));
136  }
137  std::vector<double> error(dim);
138  for(int idim=0; idim<dim; idim++)
139  error[idim] = std::abs(function_check[idim] - x_phys[idim]);
140  if (error[0] > 1e-13) {
141  std::cout << "Large error " << error[0] << std::endl;
142  for(int idim=0;idim<dim; idim++)
143  std::cout << "dim " << idim << " xref " << x_ref[idim] << " x_phys " << x_phys[idim] << " function Check " << function_check[idim] << " Error " << error[idim] << " Flag " << flag << std::endl;
144  }
145 
146  return x_ref;
147 }
148 
149 template<int dim>
150 dealii::Point<dim> CurvManifold<dim>::push_forward(const dealii::Point<dim> &chart_point) const
151 {
152  const double pi = atan(1)*4.0;
153 
154  dealii::Point<dim> x_ref;
155  dealii::Point<dim> x_phys;
156  for(int idim=0; idim<dim; idim++)
157  x_ref[idim] = chart_point[idim];
158  double beta = 1.0/10.0;
159  double alpha = 1.0/10.0;
160  if(dim==2){
161  x_phys[0] = x_ref[0] + beta*std::cos(pi/2.0*x_ref[0])*std::cos(3.0*pi/2.0*x_ref[1]);
162  x_phys[1] = x_ref[1] + beta*std::sin(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
163  }
164  else{
165  x_phys[0] =x_ref[0] + alpha*(std::cos(pi * x_ref[2]) + std::cos(pi * x_ref[1]));
166  x_phys[1] =x_ref[1] + alpha*exp(1.0-x_ref[1])*(std::sin(pi * x_ref[0]) + std::sin(pi* x_ref[2]));
167  x_phys[2] =x_ref[2] + 1.0/20.0*( std::sin(2.0 * pi * x_ref[0]) + std::sin(2.0 * pi * x_ref[1]));
168  }
169  return dealii::Point<dim> (x_phys); // Trigonometric
170 }
171 
172 template<int dim>
173 dealii::DerivativeForm<1,dim,dim> CurvManifold<dim>::push_forward_gradient(const dealii::Point<dim> &chart_point) const
174 {
175  const double pi = atan(1)*4.0;
176  dealii::DerivativeForm<1, dim, dim> dphys_dref;
177  double beta = 1.0/10.0;
178  double alpha = 1.0/10.0;
179  dealii::Point<dim> x_ref;
180  for(int idim=0; idim<dim; idim++){
181  x_ref[idim] = chart_point[idim];
182  }
183 
184  if(dim==2){
185  dphys_dref[0][0] = 1.0 - beta*pi/2.0 * std::sin(pi/2.0*x_ref[0])*std::cos(3.0*pi/2.0*x_ref[1]);
186  dphys_dref[0][1] = - beta*3.0*pi/2.0 * std::cos(pi/2.0*x_ref[0])*std::sin(3.0*pi/2.0*x_ref[1]);
187 
188  dphys_dref[1][0] = beta*2.0*pi*std::cos(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
189  dphys_dref[1][1] = 1.0 -beta*pi/2.0*std::sin(2.0*pi*(x_ref[0]))*std::sin(pi/2.0*x_ref[1]);
190  }
191  else{
192  dphys_dref[0][0] = 1.0;
193  dphys_dref[0][1] = - alpha*pi*std::sin(pi*x_ref[1]);
194  dphys_dref[0][2] = - alpha*pi*std::sin(pi*x_ref[2]);
195 
196  dphys_dref[1][0] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[0]);
197  dphys_dref[1][1] = 1.0 -alpha*exp(1.0-x_ref[1])*(std::sin(pi*x_ref[0])+std::sin(pi*x_ref[2]));
198  dphys_dref[1][2] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[2]);
199  dphys_dref[2][0] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[0]);
200  dphys_dref[2][1] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[1]);
201  dphys_dref[2][2] = 1.0;
202  }
203 
204  return dphys_dref;
205 }
206 
207 template<int dim>
208 std::unique_ptr<dealii::Manifold<dim,dim> > CurvManifold<dim>::clone() const
209 {
210  return std::make_unique<CurvManifold<dim>>();
211 }
212 
213 template <int dim>
214 static dealii::Point<dim> warp (const dealii::Point<dim> &p)
215 {
216  const double pi = atan(1)*4.0;
217  dealii::Point<dim> q = p;
218 
219  double beta =1.0/10.0;
220  double alpha =1.0/10.0;
221  if (dim == 2){
222  q[dim-2] =p[dim-2] + beta*std::cos(pi/2.0 * p[dim-2]) * std::cos(3.0 * pi/2.0 * p[dim-1]);
223  q[dim-1] =p[dim-1] + beta*std::sin(2.0 * pi * (p[dim-2])) * std::cos(pi /2.0 * p[dim-1]);
224  }
225  if(dim==3){
226  q[0] =p[0] + alpha*(std::cos(pi * p[2]) + std::cos(pi * p[1]));
227  q[1] =p[1] + alpha*exp(1.0-p[1])*(std::sin(pi * p[0]) + std::sin(pi* p[2]));
228  q[2] =p[2] + 1.0/20.0*( std::sin(2.0 * pi * p[0]) + std::sin(2.0 * pi * p[1]));
229  }
230 
231  return q;
232 }
233 
234 /****************************
235  * End of Curvilinear Grid
236  * ***************************/
237 
238 int main (int argc, char * argv[])
239 {
240 
241  dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
242  using real = double;
243  using namespace PHiLiP;
244  std::cout << std::setprecision(std::numeric_limits<long double>::digits10 + 1) << std::scientific;
245  const int dim = PHILIP_DIM;
246  const int nstate = 1;
247  dealii::ParameterHandler parameter_handler;
248  PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler);
249  dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
250 
251  PHiLiP::Parameters::AllParameters all_parameters_new;
252  all_parameters_new.parse_parameters (parameter_handler);
253 
254  // all_parameters_new.flux_nodes_type = Parameters::AllParameters::FluxNodes::GLL;
255 
256  //unsigned int poly_degree = 3;
257  double left = 0.0;
258  double right = 1.0;
259  const bool colorize = true;
260  //Generate a standard grid
261 
262  using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
263  std::shared_ptr<Triangulation> grid = std::make_shared<Triangulation>(
264  MPI_COMM_WORLD,
265  typename dealii::Triangulation<dim>::MeshSmoothing(
266  dealii::Triangulation<dim>::smoothing_on_refinement |
267  dealii::Triangulation<dim>::smoothing_on_coarsening));
268  dealii::GridGenerator::hyper_cube (*grid, left, right, colorize);
269  grid->refine_global(0);
270 
271  //Warp the grid
272  //IF WANT NON-WARPED GRID COMMENT UNTIL SAYS "NOT COMMENT"
273  dealii::GridTools::transform (&warp<dim>, *grid);
274 
275  // Assign a manifold to have curved geometry
276  const CurvManifold<dim> curv_manifold;
277  unsigned int manifold_id=0; // top face, see GridGenerator::hyper_rectangle, colorize=true
278  grid->reset_all_manifolds();
279  grid->set_all_manifold_ids(manifold_id);
280  grid->set_manifold ( manifold_id, curv_manifold );
281  //"END COMMENT" TO NOT WARP GRID
282 
283  double surf_int = 0.0;
284  for(unsigned int poly_degree = 2; poly_degree<6; poly_degree++){
285  unsigned int grid_degree = poly_degree;
286  // setup operator
287  // setup DG
288  std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(&all_parameters_new, poly_degree, poly_degree, grid_degree, grid);
289  dg->allocate_system ();
290 
291  const dealii::FESystem<dim> &fe_metric = (dg->high_order_grid->fe_system);
292  const unsigned int n_metric_dofs = fe_metric.dofs_per_cell;
293  auto metric_cell = dg->high_order_grid->dof_handler_grid.begin_active();
294 
295  dealii::QGaussLobatto<1> grid_quad(grid_degree +1);
296  const dealii::FE_DGQ<1> fe_grid(grid_degree);
297  const dealii::FESystem<1,1> fe_sys_grid(fe_grid, nstate);
298  dealii::QGauss<1> flux_quad(poly_degree +1);
299  dealii::QGauss<0> flux_quad_face(poly_degree +1);
300 
301  PHiLiP::OPERATOR::mapping_shape_functions<dim,2*dim,real> mapping_basis(nstate,poly_degree,grid_degree);
302  mapping_basis.build_1D_shape_functions_at_grid_nodes(fe_sys_grid, grid_quad);
303  mapping_basis.build_1D_shape_functions_at_flux_nodes(fe_sys_grid, flux_quad, flux_quad_face);
304 
305  PHiLiP::OPERATOR::metric_operators<real,dim,2*dim> metric_oper(nstate, poly_degree, grid_degree);
306 
307  for (auto current_cell = dg->dof_handler.begin_active(); current_cell!=dg->dof_handler.end(); ++current_cell, ++metric_cell) {
308  if (!current_cell->is_locally_owned()) continue;
309 
310  std::vector<dealii::types::global_dof_index> current_metric_dofs_indices(n_metric_dofs);
311  metric_cell->get_dof_indices (current_metric_dofs_indices);
312  std::array<std::vector<real>,dim> mapping_support_points;
313  for(int idim=0; idim<dim; idim++){
314  mapping_support_points[idim].resize(n_metric_dofs/dim);
315  }
316  dealii::QGaussLobatto<dim> vol_GLL(grid_degree +1);
317  for (unsigned int igrid_node = 0; igrid_node< n_metric_dofs/dim; ++igrid_node) {
318  for (unsigned int idof = 0; idof< n_metric_dofs; ++idof) {
319  const real val = (dg->high_order_grid->volume_nodes[current_metric_dofs_indices[idof]]);
320  const unsigned int istate = fe_metric.system_to_component_index(idof).first;
321  mapping_support_points[istate][igrid_node] += val * fe_metric.shape_value_component(idof,vol_GLL.point(igrid_node),istate);
322  }
323  }
324  const unsigned int n_quad_face_pts = dg->face_quadrature_collection[poly_degree].size();
325  const std::vector<real> &quad_weights = dg->face_quadrature_collection[poly_degree].get_weights ();
326  for (unsigned int iface=0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface) {
327 
328  metric_oper.build_facet_metric_operators(
329  iface,
330  n_quad_face_pts, n_metric_dofs/dim,
331  mapping_support_points,
332  mapping_basis,
333  false);
334 
335  const dealii::Tensor<1,dim,real> unit_normal_int = dealii::GeometryInfo<dim>::unit_normal_vector[iface];
336  std::vector<dealii::Tensor<1,dim,real>> normals_int(n_quad_face_pts);
337  for(unsigned int iquad=0; iquad<n_quad_face_pts; iquad++){
338  for(unsigned int idim=0; idim<dim; idim++){
339  normals_int[iquad][idim] = 0.0;
340  for(int idim2=0; idim2<dim; idim2++){
341  normals_int[iquad][idim] += unit_normal_int[idim2] * metric_oper.metric_cofactor_surf[idim][idim2][iquad];//\hat{n}^r * C_m^T
342  }
343  }
344  }
345  for(unsigned int iquad=0; iquad<n_quad_face_pts; iquad++){
346  for(int idim=0; idim<dim; idim++){
347  surf_int += 1.0 * quad_weights[iquad] * normals_int[iquad][idim] * 1.0;
348  }
349  }
350  }//end of face loop
351  }//end of cell loop
352  }//end poly degree loop
353  const double surf_int_mpi= (dealii::Utilities::MPI::max(surf_int, MPI_COMM_WORLD));
354 
355  if( std::abs(surf_int_mpi) > 1e-13){
356  pcout<<" Metrics Do NOT Satisfy GCL Condition\n"<<std::endl;
357  return 1;
358  }
359  else{
360  pcout<<" Metrics Satisfy GCL Condition\n"<<std::endl;
361  return 0;
362  }
363 }//end of main
CurvManifold::pull_back
virtual dealii::Point< dim > pull_back(const dealii::Point< dim > &space_point) const override
See dealii::Manifold.
Definition: GCL_Collocated_test.cpp:67
PHiLiP
Files for the baseline physics.
Definition: ADTypes.hpp:10
PHiLiP::Parameters::AllParameters
Main parameter class that contains the various other sub-parameter classes.
Definition: all_parameters.h:33
PHiLiP::DGFactory::create_discontinuous_galerkin
static std::shared_ptr< DGBase< dim, real, MeshType > > create_discontinuous_galerkin(const Parameters::AllParameters *const parameters_input, const unsigned int degree, const unsigned int max_degree_input, const unsigned int grid_degree_input, const std::shared_ptr< Triangulation > triangulation_input)
Creates a derived object DG, but returns it as DGBase.
Definition: dg_factory.cpp:10
PHiLiP::OPERATOR::metric_operators
Base metric operators class that stores functions used in both the volume and on surface.
Definition: operators.h:1086
PHiLiP::OPERATOR::mapping_shape_functions
The mapping shape functions evaluated at the desired nodes (facet set included in volume grid nodes f...
Definition: operators.h:1026
CurvManifold::clone
virtual std::unique_ptr< dealii::Manifold< dim, dim > > clone() const override
See dealii::Manifold.
Definition: GCL_Collocated_test.cpp:209
CurvManifold::push_forward_gradient
virtual dealii::DerivativeForm< 1, dim, dim > push_forward_gradient(const dealii::Point< dim > &chart_point) const override
See dealii::Manifold.
Definition: GCL_Collocated_test.cpp:174
PHiLiP::Parameters::AllParameters::parse_parameters
void parse_parameters(dealii::ParameterHandler &prm)
Retrieve parameters from dealii::ParameterHandler.
Definition: all_parameters.cpp:377
CurvManifold::push_forward
virtual dealii::Point< dim > push_forward(const dealii::Point< dim > &chart_point) const override
See dealii::Manifold.
Definition: GCL_Collocated_test.cpp:151
PHiLiP::Parameters::AllParameters::declare_parameters
static void declare_parameters(dealii::ParameterHandler &prm)
Declare parameters that can be set as inputs and set up the default options.
Definition: all_parameters.cpp:34