[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
weight_adjusted_mass_inverse_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 //#include <ctime>
9 #include <time.h>
10 
11 #include <deal.II/distributed/solution_transfer.h>
12 
13 #include "testing/tests.h"
14 
15 #include<fstream>
16 #include <deal.II/base/parameter_handler.h>
17 #include <deal.II/base/tensor.h>
18 #include <deal.II/numerics/vector_tools.h>
19 
20 #include <deal.II/grid/grid_generator.h>
21 #include <deal.II/grid/grid_refinement.h>
22 #include <deal.II/grid/grid_tools.h>
23 #include <deal.II/grid/grid_out.h>
24 #include <deal.II/grid/grid_in.h>
25 
26 #include <deal.II/dofs/dof_handler.h>
27 #include <deal.II/dofs/dof_tools.h>
28 #include <deal.II/dofs/dof_renumbering.h>
29 
30 #include <deal.II/dofs/dof_accessor.h>
31 
32 #include <deal.II/lac/vector.h>
33 #include <deal.II/lac/dynamic_sparsity_pattern.h>
34 #include <deal.II/lac/sparse_matrix.h>
35 
36 #include <deal.II/meshworker/dof_info.h>
37 
38 #include <deal.II/base/convergence_table.h>
39 
40 // Finally, we take our exact solution from the library as well as volume_quadrature
41 // and additional tools.
42 #include <deal.II/numerics/data_out.h>
43 #include <deal.II/numerics/data_out_dof_data.h>
44 #include <deal.II/numerics/vector_tools.h>
45 #include <deal.II/numerics/vector_tools.templates.h>
46 
47 #include "parameters/all_parameters.h"
48 #include "parameters/parameters.h"
49 #include "dg/dg_base.hpp"
50 #include <deal.II/grid/manifold_lib.h>
51 #include <deal.II/fe/mapping_q.h>
52 #include "dg/dg_factory.hpp"
53 #include "operators/operators.h"
54 //#include <GCL_test.h>
55 
56 const double TOLERANCE = 1E-6;
57 using namespace std;
58 //namespace PHiLiP {
59 
60 template <int dim>
61 class CurvManifold: public dealii::ChartManifold<dim,dim,dim> {
62  virtual dealii::Point<dim> pull_back(const dealii::Point<dim> &space_point) const override;
63  virtual dealii::Point<dim> push_forward(const dealii::Point<dim> &chart_point) const override;
64  virtual dealii::DerivativeForm<1,dim,dim> push_forward_gradient(const dealii::Point<dim> &chart_point) const override;
65 
66  virtual std::unique_ptr<dealii::Manifold<dim,dim> > clone() const override;
67 };
68 
69 template<int dim>
70 dealii::Point<dim> CurvManifold<dim>::pull_back(const dealii::Point<dim> &space_point) const
71 {
72  using namespace PHiLiP;
73  const double pi = atan(1)*4.0;
74  dealii::Point<dim> x_ref;
75  dealii::Point<dim> x_phys;
76  for(int idim=0; idim<dim; idim++){
77  x_ref[idim] = space_point[idim];
78  x_phys[idim] = space_point[idim];
79  }
80  dealii::Vector<double> function(dim);
81  dealii::FullMatrix<double> derivative(dim);
82  double beta =1.0/10.0;
83  double alpha =1.0/10.0;
84  int flag =0;
85  while(flag != dim){
86  if(dim==2){
87  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]);
88  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]);
89  }
90  else{
91  function[0] = x_ref[0] - x_phys[0] +alpha*(std::cos(pi * x_ref[2]) + std::cos(pi * x_ref[1]));
92  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]));
93  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]));
94  }
95 
96 
97  if(dim==2){
98  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]);
99  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]);
100 
101  derivative[1][0] = beta*2.0*pi*std::cos(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
102  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]);
103  }
104  else{
105  derivative[0][0] = 1.0;
106  derivative[0][1] = - alpha*pi*std::sin(pi*x_ref[1]);
107  derivative[0][2] = - alpha*pi*std::sin(pi*x_ref[2]);
108 
109  derivative[1][0] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[0]);
110  derivative[1][1] = 1.0 -alpha*exp(1.0-x_ref[1])*(std::sin(pi*x_ref[0])+std::sin(pi*x_ref[2]));
111  derivative[1][2] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[2]);
112  derivative[2][0] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[0]);
113  derivative[2][1] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[1]);
114  derivative[2][2] = 1.0;
115  }
116 
117  dealii::FullMatrix<double> Jacobian_inv(dim);
118  Jacobian_inv.invert(derivative);
119  dealii::Vector<double> Newton_Step(dim);
120  Jacobian_inv.vmult(Newton_Step, function);
121  for(int idim=0; idim<dim; idim++){
122  x_ref[idim] -= Newton_Step[idim];
123  }
124  flag=0;
125  for(int idim=0; idim<dim; idim++){
126  if(std::abs(function[idim]) < 1e-15)
127  flag++;
128  }
129  if(flag == dim)
130  break;
131  }
132  std::vector<double> function_check(dim);
133  if(dim==2){
134  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]);
135  function_check[1] = x_ref[1] + beta*std::sin(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
136  }
137  else{
138  function_check[0] = x_ref[0] +alpha*(std::cos(pi * x_ref[2]) + std::cos(pi * x_ref[1]));
139  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]));
140  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]));
141  }
142  std::vector<double> error(dim);
143  for(int idim=0; idim<dim; idim++)
144  error[idim] = std::abs(function_check[idim] - x_phys[idim]);
145  if (error[0] > 1e-13) {
146  std::cout << "Large error " << error[0] << std::endl;
147  for(int idim=0;idim<dim; idim++)
148  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;
149  }
150 
151  return x_ref;
152 }
153 
154 template<int dim>
155 dealii::Point<dim> CurvManifold<dim>::push_forward(const dealii::Point<dim> &chart_point) const
156 {
157  const double pi = atan(1)*4.0;
158 
159  dealii::Point<dim> x_ref;
160  dealii::Point<dim> x_phys;
161  for(int idim=0; idim<dim; idim++)
162  x_ref[idim] = chart_point[idim];
163  double beta = 1.0/10.0;
164  double alpha = 1.0/10.0;
165  if(dim==2){
166  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]);
167  x_phys[1] = x_ref[1] + beta*std::sin(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
168  }
169  else{
170  x_phys[0] =x_ref[0] + alpha*(std::cos(pi * x_ref[2]) + std::cos(pi * x_ref[1]));
171  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]));
172  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]));
173  }
174  return dealii::Point<dim> (x_phys); // Trigonometric
175 }
176 
177 template<int dim>
178 dealii::DerivativeForm<1,dim,dim> CurvManifold<dim>::push_forward_gradient(const dealii::Point<dim> &chart_point) const
179 {
180  const double pi = atan(1)*4.0;
181  dealii::DerivativeForm<1, dim, dim> dphys_dref;
182  double beta = 1.0/10.0;
183  double alpha = 1.0/10.0;
184  dealii::Point<dim> x_ref;
185  for(int idim=0; idim<dim; idim++){
186  x_ref[idim] = chart_point[idim];
187  }
188 
189  if(dim==2){
190  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]);
191  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]);
192 
193  dphys_dref[1][0] = beta*2.0*pi*std::cos(2.0*pi*(x_ref[0]))*std::cos(pi/2.0*x_ref[1]);
194  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]);
195  }
196  else{
197  dphys_dref[0][0] = 1.0;
198  dphys_dref[0][1] = - alpha*pi*std::sin(pi*x_ref[1]);
199  dphys_dref[0][2] = - alpha*pi*std::sin(pi*x_ref[2]);
200 
201  dphys_dref[1][0] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[0]);
202  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]));
203  dphys_dref[1][2] = alpha*pi*exp(1.0-x_ref[1])*std::cos(pi*x_ref[2]);
204  dphys_dref[2][0] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[0]);
205  dphys_dref[2][1] = 1.0/10.0*pi*std::cos(2.0*pi*x_ref[1]);
206  dphys_dref[2][2] = 1.0;
207  }
208 
209  return dphys_dref;
210 }
211 
212 template<int dim>
213 std::unique_ptr<dealii::Manifold<dim,dim> > CurvManifold<dim>::clone() const
214 {
215  return std::make_unique<CurvManifold<dim>>();
216 }
217 
218 template <int dim>
219 static dealii::Point<dim> warp (const dealii::Point<dim> &p)
220 {
221  const double pi = atan(1)*4.0;
222  dealii::Point<dim> q = p;
223 
224  double beta =1.0/10.0;
225  double alpha =1.0/10.0;
226  if (dim == 2){
227  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]);
228  q[dim-1] =p[dim-1] + beta*std::sin(2.0 * pi * (p[dim-2])) * std::cos(pi /2.0 * p[dim-1]);
229  }
230  if(dim==3){
231  q[0] =p[0] + alpha*(std::cos(pi * p[2]) + std::cos(pi * p[1]));
232  q[1] =p[1] + alpha*exp(1.0-p[1])*(std::sin(pi * p[0]) + std::sin(pi* p[2]));
233  q[2] =p[2] + 1.0/20.0*( std::sin(2.0 * pi * p[0]) + std::sin(2.0 * pi * p[1]));
234  }
235 
236  return q;
237 }
238 
239 /****************************
240  * End of Curvilinear Grid
241  * ***************************/
242 
243 template <int dim, typename real>
244 void compute_inverse_mass_matrix(
245  std::shared_ptr < PHiLiP::OPERATOR::OperatorsBase<dim,real,2*dim> > &operators,
246  const std::array<std::vector<real>,PHILIP_DIM> &mapping_support_points,
247  const unsigned int n_metric_dofs,
248  const unsigned int n_quad_pts, const unsigned int n_dofs_cell,
249  const unsigned int poly_degree, const unsigned int grid_degree,
250  const std::vector<real> &quad_weights,
251  dealii::FullMatrix<real> &mass_inv)
252 {
253  std::vector<real> determinant_Jacobian(n_quad_pts);
254  operators->build_local_vol_determinant_Jac(grid_degree, poly_degree, n_quad_pts, n_metric_dofs, mapping_support_points, determinant_Jacobian);
255 
256  std::vector<real> JxW(n_quad_pts);
257  for(unsigned int iquad=0; iquad<n_quad_pts; iquad++){
258  JxW[iquad] = quad_weights[iquad] * determinant_Jacobian[iquad];
259  }
260  dealii::FullMatrix<real> local_mass_matrix(n_dofs_cell);
261  operators->build_local_Mass_Matrix(JxW, n_dofs_cell, n_quad_pts, poly_degree, local_mass_matrix);
262 
263  //For flux reconstruction
264  dealii::FullMatrix<real> Flux_Reconstruction_operator(n_dofs_cell);
265  operators->build_local_Flux_Reconstruction_operator(local_mass_matrix, n_dofs_cell, poly_degree, Flux_Reconstruction_operator);
266  for (unsigned int itest=0; itest<n_dofs_cell; ++itest) {
267  for (unsigned int itrial=0; itrial<n_dofs_cell; ++itrial) {
268  local_mass_matrix[itest][itrial] = local_mass_matrix[itest][itrial] + Flux_Reconstruction_operator[itest][itrial];
269  }
270  }
271 
272  mass_inv.invert(local_mass_matrix);
273 }
274 
275 template <int dim, typename real>
276 void compute_weighted_inverse_mass_matrix(std::shared_ptr < PHiLiP::OPERATOR::OperatorsBase<dim,real,2*dim> > &operators,
277  const std::array<std::vector<real>,PHILIP_DIM> &mapping_support_points, const unsigned int n_metric_dofs,
278  const unsigned int n_quad_pts, const unsigned int n_dofs_cell,
279  const unsigned int poly_degree, const unsigned int grid_degree,
280  const std::vector<real> &quad_weights,
281  dealii::FullMatrix<real> &mass_inv)
282 {
283  std::vector<real> determinant_Jacobian(n_quad_pts);
284  operators->build_local_vol_determinant_Jac(grid_degree, poly_degree, n_quad_pts, n_metric_dofs, mapping_support_points, determinant_Jacobian);
285 
286  std::vector<real> W_J_inv(n_quad_pts);
287  for(unsigned int iquad=0; iquad<n_quad_pts; iquad++){
288  W_J_inv[iquad] = quad_weights[iquad] / determinant_Jacobian[iquad];
289  }
290  dealii::FullMatrix<real> local_mass_matrix(n_dofs_cell);
291  operators->build_local_Mass_Matrix(W_J_inv, n_dofs_cell, n_quad_pts, poly_degree, local_mass_matrix);
292  // For flux reconstruction
293  dealii::FullMatrix<real> Flux_Reconstruction_operator(n_dofs_cell);
294  operators->build_local_Flux_Reconstruction_operator(local_mass_matrix, n_dofs_cell, poly_degree, Flux_Reconstruction_operator);
295  for (unsigned int itest=0; itest<n_dofs_cell; ++itest) {
296  for (unsigned int itrial=0; itrial<n_dofs_cell; ++itrial) {
297  // local_mass_matrix[itest][itrial] = local_mass_matrix[itest][itrial] + Flux_Reconstruction_operator[itest][itrial];
298  mass_inv[itest][itrial] = local_mass_matrix[itest][itrial] + Flux_Reconstruction_operator[itest][itrial];
299  }
300  }
301 }
302 
303 /*******************************
304  * END OF MASS INV FUNCTIONS
305  * ****************************/
306 
307 int main (int argc, char * argv[])
308 {
309  dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
310  using real = double;
311  using namespace PHiLiP;
312  std::cout << std::setprecision(std::numeric_limits<long double>::digits10 + 1) << std::scientific;
313  const int dim = PHILIP_DIM;
314  dealii::ParameterHandler parameter_handler;
315  PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler);
316  dealii::ConditionalOStream pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0);
317 
318  PHiLiP::Parameters::AllParameters all_parameters_new;
319  all_parameters_new.parse_parameters (parameter_handler);
320 
321  // all_parameters_new.flux_nodes_type = Parameters::AllParameters::FluxNodes::GLL;
322  all_parameters_new.use_curvilinear_split_form=true;
323  all_parameters_new.flux_reconstruction_type = Parameters::AllParameters::Flux_Reconstruction::cPlus;
324 
325  // unsigned int poly_degree = 3;
326  double left = 0.0;
327  double right = 1.0;
328  const bool colorize = true;
329  dealii::ConvergenceTable convergence_table;
330  const unsigned int igrid_start = 0;
331  const unsigned int n_grids = 1;
332  // setup time
333  // time_t tstart=0, tend=0, tstart_weight=0, tend_weight=0;
334  clock_t time_normal, time_weighted;
335 
336  for(unsigned int poly_degree = 6; poly_degree<7; poly_degree++){
337  unsigned int grid_degree = poly_degree;
338  for(unsigned int igrid=igrid_start; igrid<n_grids; ++igrid){
339  pcout<<" Grid Index"<<igrid<<std::endl;
340 
341  //Generate a standard grid
342  using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
343  std::shared_ptr<Triangulation> grid = std::make_shared<Triangulation>(
344  MPI_COMM_WORLD,
345  typename dealii::Triangulation<dim>::MeshSmoothing(
346  dealii::Triangulation<dim>::smoothing_on_refinement |
347  dealii::Triangulation<dim>::smoothing_on_coarsening));
348  dealii::GridGenerator::hyper_cube (*grid, left, right, colorize);
349  grid->refine_global(igrid);
350  pcout<<" made grid for Index"<<igrid<<std::endl;
351 
352  //Warp the grid
353  //IF WANT NON-WARPED GRID COMMENT UNTIL SAYS "NOT COMMENT"
354  dealii::GridTools::transform (&warp<dim>, *grid);
355 
356  // Assign a manifold to have curved geometry
357  const CurvManifold<dim> curv_manifold;
358  unsigned int manifold_id=0; // top face, see GridGenerator::hyper_rectangle, colorize=true
359  grid->reset_all_manifolds();
360  grid->set_all_manifold_ids(manifold_id);
361  grid->set_manifold ( manifold_id, curv_manifold );
362  //"END COMMENT" TO NOT WARP GRID
363 
364  // setup operator
365  // OPERATOR::OperatorsBase<dim,real> operators(&all_parameters_new, nstate, poly_degree, poly_degree, grid_degree);
366  // OPERATOR::OperatorsBaseState<dim,real,nstate,2*dim> operators(&all_parameters_new, poly_degree, poly_degree);
367  // setup DG
368  // std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(&all_parameters_new, poly_degree, grid);
369  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);
370  dg->allocate_system ();
371 
372  dealii::IndexSet locally_owned_dofs;
373  dealii::IndexSet ghost_dofs;
374  dealii::IndexSet locally_relevant_dofs;
375  locally_owned_dofs = dg->dof_handler.locally_owned_dofs();
376  dealii::DoFTools::extract_locally_relevant_dofs(dg->dof_handler, ghost_dofs);
377  locally_relevant_dofs = ghost_dofs;
378  ghost_dofs.subtract_set(locally_owned_dofs);
379 
380  // setup metric and solve
381  const unsigned int max_dofs_per_cell = dg->dof_handler.get_fe_collection().max_dofs_per_cell();
382  std::vector<dealii::types::global_dof_index> current_dofs_indices(max_dofs_per_cell);
383  const unsigned int n_dofs_cell = dg->operators->fe_collection_basis[poly_degree].dofs_per_cell;
384  const unsigned int n_quad_pts = dg->operators->volume_quadrature_collection[poly_degree].size();
385 
386  const dealii::FESystem<dim> &fe_metric = (dg->high_order_grid->fe_system);
387  const unsigned int n_metric_dofs = fe_metric.dofs_per_cell;
388  auto metric_cell = dg->high_order_grid->dof_handler_grid.begin_active();
389 
390  //loop over cells and do normal inv
391  pcout<<"time to do normal"<<std::endl;
392  // tstart = time(0);
393  time_normal = clock();
394  for (auto current_cell = dg->dof_handler.begin_active(); current_cell!=dg->dof_handler.end(); ++current_cell, ++metric_cell) {
395  if (!current_cell->is_locally_owned()) continue;
396 
397  // pcout<<"grid degree "<<grid_degree<<" metric dofs "<<n_metric_dofs<<std::endl;
398  std::vector<dealii::types::global_dof_index> current_metric_dofs_indices(n_metric_dofs);
399  metric_cell->get_dof_indices (current_metric_dofs_indices);
400  std::array<std::vector<real>,dim> mapping_support_points;
401  for(int idim=0; idim<dim; idim++){
402  mapping_support_points[idim].resize(n_metric_dofs/dim);
403  }
404  dealii::QGaussLobatto<dim> vol_GLL(grid_degree +1);
405  for (unsigned int igrid_node = 0; igrid_node< n_metric_dofs/dim; ++igrid_node) {
406  for (unsigned int idof = 0; idof< n_metric_dofs; ++idof) {
407  const real val = (dg->high_order_grid->volume_nodes[current_metric_dofs_indices[idof]]);
408  const unsigned int istate = fe_metric.system_to_component_index(idof).first;
409  mapping_support_points[istate][igrid_node] += val * fe_metric.shape_value_component(idof,vol_GLL.point(igrid_node),istate);
410  }
411  }
412  const std::vector<real> &quad_weights = dg->operators->volume_quadrature_collection[poly_degree].get_weights();
413 
414  //build ESFR mass matrix and invert regularly
415  dealii::FullMatrix<real> mass_inv(n_dofs_cell);
416  time_normal = clock();
417  compute_inverse_mass_matrix(dg->operators, mapping_support_points, n_metric_dofs/dim, n_quad_pts, n_dofs_cell, poly_degree, grid_degree, quad_weights, mass_inv);
418  time_normal = clock()-time_normal;
419 
420  }//end of cell loop
421 
422  // tend = time(0);
423  // time_normal = clock()-time_normal;
424 
425  pcout<<"time to do weighted"<<std::endl;
426  metric_cell = dg->high_order_grid->dof_handler_grid.begin_active();
427  //loop over cells and do weight inv
428  // tstart_weight = time(0);
429  time_weighted = clock();
430  for (auto current_cell = dg->dof_handler.begin_active(); current_cell!=dg->dof_handler.end(); ++current_cell, ++metric_cell) {
431  if (!current_cell->is_locally_owned()) continue;
432 
433  // pcout<<"grid degree "<<grid_degree<<" metric dofs "<<n_metric_dofs<<std::endl;
434  std::vector<dealii::types::global_dof_index> current_metric_dofs_indices(n_metric_dofs);
435  metric_cell->get_dof_indices (current_metric_dofs_indices);
436  std::array<std::vector<real>,dim> mapping_support_points;
437  for(int idim=0; idim<dim; idim++){
438  mapping_support_points[idim].resize(n_metric_dofs/dim);
439  }
440  dealii::QGaussLobatto<dim> vol_GLL(grid_degree +1);
441  for (unsigned int igrid_node = 0; igrid_node< n_metric_dofs/dim; ++igrid_node) {
442  for (unsigned int idof = 0; idof< n_metric_dofs; ++idof) {
443  const real val = (dg->high_order_grid->volume_nodes[current_metric_dofs_indices[idof]]);
444  const unsigned int istate = fe_metric.system_to_component_index(idof).first;
445  mapping_support_points[istate][igrid_node] += val * fe_metric.shape_value_component(idof,vol_GLL.point(igrid_node),istate);
446  }
447  }
448  const std::vector<real> &quad_weights = dg->operators->volume_quadrature_collection[poly_degree].get_weights();
449 
450  //do weight-adjusted inverse ESFR mass matrix
451  dealii::FullMatrix<real> mass_inv(n_dofs_cell);
452  time_weighted = clock();
453  compute_weighted_inverse_mass_matrix(dg->operators, mapping_support_points, n_metric_dofs/dim, n_quad_pts, n_dofs_cell, poly_degree, grid_degree, quad_weights, mass_inv);
454  time_weighted = clock() - time_weighted;
455 
456  }//end of cell loop
457 
458  // tend_weight = time(0);
459  // time_weighted = clock() - time_weighted;
460  }//end grid refinement loop
461  }//end poly degree loop
462 
463  // pcout<<"Normal Mass inv took "<<difftime(tend, tstart)<<" seconds (s)."<<std::endl;
464  // pcout<<"Weighted Mass inv took "<<difftime(tend_weight, tstart_weight)<<" seconds (s)."<<std::endl;
465  // pcout<<"Normal Mass inv took "<<time_normal ((float)time_normal)/CLOCKS_PER_SEC<<" seconds (s)."<<std::endl;
466  printf(" it took %g seconds normal\n",((float)time_normal)/CLOCKS_PER_SEC);
467  printf(" it took %g seconds weighted\n",((float)time_weighted)/CLOCKS_PER_SEC);
468  // pcout<<"Weighted Mass inv took "<<time_weighted<<" seconds (s)."<<std::endl;
469 
470  // if(difftime(tend, tstart) < difftime(tend_weight, tstart_weight)){
471  if(time_normal < time_weighted){
472  pcout<<"Weighted inv not faster!"<<std::endl;
473  return 1;
474  }
475  else {
476  return 0;
477  }
478 }//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::Parameters::AllParameters::flux_reconstruction_type
Flux_Reconstruction flux_reconstruction_type
Store flux reconstruction type.
Definition: all_parameters.h:280
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
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
PHiLiP::OPERATOR::OperatorsBase
Operator base class.
Definition: operators.h:53
PHiLiP::Parameters::AllParameters::use_curvilinear_split_form
bool use_curvilinear_split_form
Flag to use curvilinear metric split form.
Definition: all_parameters.h:117