Flow123d  master-f44eb46
assembly_flow_output.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file assembly_flow_output.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_FLOW_OUTPUT_HH_
19 #define ASSEMBLY_FLOW_OUTPUT_HH_
20 
21 #include <fstream>
22 
26 #include "flow/darcy_flow_lmh.hh"
28 #include "fields/field_fe.hh"
29 
30 #include "mesh/mesh.h"
31 #include "mesh/accessors.hh"
32 #include "fem/fe_p.hh"
33 #include "fem/fe_values.hh"
34 #include "fem/fe_rt.hh"
35 #include "fem/fe_values_views.hh"
36 #include "fem/fe_system.hh"
38 
39 
40 /**
41  * Calculate approximation of L2 norm for:
42  * 1) difference between regularized pressure and analytical solution (using FunctionPython)
43  * 2) difference between RT velocities and analytical solution
44  * 3) difference of divergence
45  *
46  * TODO:
47  * 1) implement field objects
48  * 2) implement DG_P2 finite elements
49  * 3) implement pressure postprocessing (result is DG_P2 field)
50  * 4) implement calculation of L2 norm for two field (compute the norm and values on individual elements as P0 field)
51  *
52  */
53 template <unsigned int dim>
55 {
56 public:
57  typedef typename DarcyLMH::EqFields EqFields;
59 
60  static constexpr const char * name() { return "L2DifferenceAssembly"; }
61 
62  /// Constructor.
63  L2DifferenceAssembly(EqFields *eq_fields, EqData *eq_data)
64  : AssemblyBase<dim>(2), eq_fields_(eq_fields), eq_data_(eq_data) {
71  }
72 
73  /// Destructor.
74  virtual ~L2DifferenceAssembly() {}
75 
76  /// Initialize auxiliary vectors and other data members
77  void initialize(ElementCacheMap *element_cache_map) {
78  this->element_cache_map_ = element_cache_map;
79 
81  fe_p0_ = std::make_shared< FE_P_disc<dim> >(0);
82  fe_values_.initialize(*this->quad_, *fe_p0_, flags);
83  fv_rt_.initialize(*this->quad_, fe_rt_, flags);
84 
85  fluxes_.resize(dim+1);
86  }
87 
88 
89  /// Assemble integral over element
90  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
91  {
92  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
93 
94  ElementAccessor<3> ele = cell.elm();
95  fv_rt_.reinit(ele);
96  fe_values_.reinit(ele);
97 
98  // get coefficients on the current element
99  for (unsigned int li = 0; li < ele->n_sides(); li++) {
100  fluxes_[li] = eq_data_->flow_data_->full_solution.get( cell.get_loc_dof_indices()[li] );
101  //pressure_traces[li] = diff_eq_data_->dh->side_scalar( *(ele->side( li ) ) );
102  }
103  // TODO: replace with DHCell getter when available for FESystem component
104  double pressure_mean = eq_data_->flow_data_->full_solution.get( cell.get_loc_dof_indices()[ cell.n_dofs()/2 ] );
105 
106  double velocity_diff=0, divergence_diff=0, pressure_diff=0, diff;
107 
108  // 1d: mean_x_squared = 1/6 (v0^2 + v1^2 + v0*v1)
109  // 2d: mean_x_squared = 1/12 (v0^2 + v1^2 +v2^2 + v0*v1 + v0*v2 + v1*v2)
110  double mean_x_squared=0;
111  for(unsigned int i_node=0; i_node < ele->n_nodes(); i_node++ )
112  for(unsigned int j_node=0; j_node < ele->n_nodes(); j_node++ )
113  {
114  mean_x_squared += (i_node == j_node ? 2.0 : 1.0) / ( 6 * dim ) // multiply by 2 on diagonal
115  * arma::dot( *ele.node(i_node), *ele.node(j_node));
116  }
117 
118  unsigned int i_point=0, oposite_node;
119  for (auto p : this->bulk_points(element_patch_idx) )
120  {
121  q_point_ = fe_values_.point(i_point);
122 
128 
129  // compute postprocesed pressure
130  diff = 0;
131  for(unsigned int i_shape=0; i_shape < ele->n_sides(); i_shape++) {
132  oposite_node = RefElement<dim>::oposite_node(i_shape);
133 
134  diff += fluxes_[ i_shape ] *
135  ( arma::dot( q_point_, q_point_ )/ 2
136  - mean_x_squared / 2
137  - arma::dot( q_point_, *ele.node(oposite_node) )
138  + arma::dot( ele.centre(), *ele.node(oposite_node) )
139  );
140  }
141 
142  diff = - (1.0 / conductivity_) * diff / dim / ele.measure() / cross_ + pressure_mean ;
143  diff = ( diff - ref_pressure_);
144  pressure_diff += diff * diff * fe_values_.JxW(i_point);
145 
146  // velocity difference
147  flux_in_q_point_.zeros();
148  for(unsigned int i_shape=0; i_shape < ele->n_sides(); i_shape++) {
149  flux_in_q_point_ += fluxes_[ i_shape ]
150  * fv_rt_.vector_view(0).value(i_shape, i_point)
151  / cross_;
152  }
153 
155  velocity_diff += dot(flux_in_q_point_, flux_in_q_point_) * fe_values_.JxW(i_point);
156 
157  // divergence diff
158  diff = 0;
159  for(unsigned int i_shape=0; i_shape < ele->n_sides(); i_shape++) diff += fluxes_[ i_shape ];
160  diff = ( diff / ele.measure() / cross_ - ref_divergence_);
161  divergence_diff += diff * diff * fe_values_.JxW(i_point);
162 
163  i_point++;
164  }
165 
166  // DHCell constructed with diff fields DH, get DOF indices of actual element
167  DHCellAccessor sub_dh_cell = cell.cell_with_other_dh(eq_data_->dh_.get());
168  IntIdx idx = sub_dh_cell.get_loc_dof_indices()[0];
169 
170  auto velocity_data = eq_data_->vel_diff_ptr->vec();
171  velocity_data.set( idx, sqrt(velocity_diff) );
172  eq_data_->velocity_error[dim-1] += velocity_diff;
173  if (dim == 2 && eq_data_->velocity_mask.size() != 0 ) {
174  eq_data_->mask_vel_error += (eq_data_->velocity_mask[ ele.idx() ])? 0 : velocity_diff;
175  }
176 
177  auto pressure_data = eq_data_->pressure_diff_ptr->vec();
178  pressure_data.set( idx, sqrt(pressure_diff) );
179  eq_data_->pressure_error[dim-1] += pressure_diff;
180 
181  auto div_data = eq_data_->div_diff_ptr->vec();
182  div_data.set( idx, sqrt(divergence_diff) );
183  eq_data_->div_error[dim-1] += divergence_diff;
184  }
185 
186  /// Implements @p AssemblyBase::begin.
187  void begin() override
188  {
189  DebugOut() << "l2 norm output\n";
190 
192  for(unsigned int j=0; j<3; j++){
193  eq_data_->pressure_error[j] = 0;
194  eq_data_->velocity_error[j] = 0;
195  eq_data_->div_error[j] = 0;
196  }
197  }
198 
199  /// Implements @p AssemblyBase::end.
200  void end() override
201  {
202  ofstream os( FilePath("solution_error", FilePath::output_file) );
203 
204  // square root for L2 norm
205  for(unsigned int j=0; j<3; j++){
208  eq_data_->div_error[j] = sqrt(eq_data_->div_error[j]);
209  }
211 
212  os << "l2 norm output\n\n"
213  << "pressure error 1d: " << eq_data_->pressure_error[0] << endl
214  << "pressure error 2d: " << eq_data_->pressure_error[1] << endl
215  << "pressure error 3d: " << eq_data_->pressure_error[2] << endl
216  << "velocity error 1d: " << eq_data_->velocity_error[0] << endl
217  << "velocity error 2d: " << eq_data_->velocity_error[1] << endl
218  << "velocity error 3d: " << eq_data_->velocity_error[2] << endl
219  << "masked velocity error 2d: " << eq_data_->mask_vel_error <<endl
220  << "div error 1d: " << eq_data_->div_error[0] << endl
221  << "div error 2d: " << eq_data_->div_error[1] << endl
222  << "div error 3d: " << eq_data_->div_error[2];
223  }
224 
225 
226 protected:
227  /// Sub field set contains fields used in calculation.
229 
230  /// Data objects shared with Flow equation
233 
234  /// following is used for calculation of postprocessed pressure difference
235  /// and comparison to analytical solution
236  std::shared_ptr<FE_P_disc<dim>> fe_p0_;
238 
239  /// FEValues for velocity.
242 
243  std::vector<double> fluxes_; ///< Precomputed fluxes on element sides.
244  arma::vec3 flux_in_q_point_; ///< Precomputed flux in quadrature point.
245  arma::vec3 q_point_; ///< Local coords of quadrature point.
246  arma::vec3 ref_flux_; ///< Field result.
247  double ref_pressure_, ref_divergence_, conductivity_, cross_; ///< Field results.
248 
249  template < template<IntDim...> class DimAssembly>
250  friend class GenericAssembly;
251 
252 };
253 
254 
255 /**
256  * Compute output of internal flow data.
257  */
258 template <unsigned int dim>
260 {
261 public:
262  typedef typename DarcyLMH::EqFields EqFields;
264 
265  static constexpr const char * name() { return "OutputInternalFlowAssembly"; }
266 
267  /// Constructor.
269  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
272  }
273 
274  /// Destructor.
276 
277  /// Initialize auxiliary vectors and other data members
278  void initialize(ElementCacheMap *element_cache_map) {
279  this->element_cache_map_ = element_cache_map;
280  }
281 
282 
283  /// Assemble integral over element
284  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
285  {
286  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
287 
288  ElementAccessor<3> ele = cell.elm();
289  LocDofVec indices = cell.get_loc_dof_indices();
290 
291  std::stringstream ss;
292  // pressure
293  ss << fmt::format("{} {} ", cell.elm().input_id(), eq_data_->flow_data_->full_solution.get(indices[ele->n_sides()]));
294 
295  // velocity at element center
296  auto p = *( this->bulk_points(element_patch_idx).begin() );
298  for (unsigned int i = 0; i < 3; i++)
299  ss << flux_in_center_[i] << " ";
300 
301  // number of sides
302  ss << ele->n_sides() << " ";
303 
304  // use node permutation to permute sides
305  auto &new_to_old_node = ele.orig_nodes_order();
306  std::vector<uint> old_to_new_side(ele->n_sides());
307  for (unsigned int i = 0; i < ele->n_sides(); i++) {
308  // According to RefElement<dim>::opposite_node()
309  uint new_opp_node = ele->n_sides() - i - 1;
310  uint old_opp_node = new_to_old_node[new_opp_node];
311  uint old_iside = ele->n_sides() - old_opp_node - 1;
312  old_to_new_side[old_iside] = i;
313  }
314 
315  // pressure on edges
316  // unsigned int lid = ele->n_sides() + 1;
317  for (unsigned int i = 0; i < ele->n_sides(); i++) {
318  uint new_lid = ele->n_sides() + 1 + old_to_new_side[i];
319  ss << eq_data_->flow_data_->full_solution.get(indices[new_lid]) << " ";
320  }
321  // fluxes on sides
322  for (unsigned int i = 0; i < ele->n_sides(); i++) {
323  uint new_iside = old_to_new_side[i];
324  ss << eq_data_->flow_data_->full_solution.get(indices[new_iside]) << " ";
325  }
326 
327  // remove last white space
328  string line = ss.str();
329  eq_data_->raw_output_strings_[cell.elm_idx()] = line.substr(0, line.size()-1);
330  }
331 
332  /// Implements @p AssemblyBase::begin.
333  void begin() override
334  {
335  DebugOut() << "Internal flow data output\n";
336 
337  eq_data_->raw_output_strings_.resize( eq_data_->flow_data_->dh_->n_own_cells() );
338  }
339 
340  /// Implements @p AssemblyBase::end.
341  void end() override
342  {
343  //char dbl_fmt[ 16 ]= "%.8g ";
344  // header
345  eq_data_->raw_output_file << "// fields:\n//ele_id ele_presure flux_in_barycenter[3] n_sides side_pressures[n] side_fluxes[n]\n";
346  eq_data_->raw_output_file << fmt::format("$FlowField\nT={}\n", eq_data_->time_->t());
347  eq_data_->raw_output_file << fmt::format("{}\n" , eq_data_->flow_data_->dh_->mesh()->n_elements() );
348 
349  auto permutation_vec = eq_data_->flow_data_->dh_->mesh()->element_permutations();
350  for (unsigned int i_elem=0; i_elem<eq_data_->flow_data_->dh_->n_own_cells(); ++i_elem) {
351  eq_data_->raw_output_file << eq_data_->raw_output_strings_[ permutation_vec[i_elem] ] << endl;
352  }
353 
354  eq_data_->raw_output_file << "$EndFlowField\n" << endl;
355  }
356 
357 
358 protected:
359  /// Sub field set contains fields used in calculation.
361 
362  /// Data objects shared with Flow equation
365 
367 
368  template < template<IntDim...> class DimAssembly>
369  friend class GenericAssembly;
370 
371 };
372 
373 
374 #endif /* ASSEMBLY_FLOW_OUTPUT_HH_ */
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
Quadrature * quad_
Quadrature used in assembling methods.
int active_integrals_
Holds mask of active integrals.
ElementCacheMap * element_cache_map_
ElementCacheMap shared with GenericAssembly object.
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
unsigned int n_dofs() const
Return number of dofs on given cell.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
DHCellAccessor cell_with_other_dh(const DOFHandlerMultiDim *dh) const
Create new accessor with same local idx and given DOF handler. Actual and given DOF handler must be c...
unsigned int dim() const
Return dimension of element appropriate to cell.
unsigned int elm_idx() const
Return serial idx to element of loc_ele_idx_.
Output specific field stuff.
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > div_diff_ptr
std::shared_ptr< SubDOFHandlerMultiDim > dh_
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > vel_diff_ptr
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > pressure_diff_ptr
std::shared_ptr< DarcyLMH::EqData > flow_data_
std::shared_ptr< DarcyLMH::EqData > flow_data_
ofstream raw_output_file
Raw data output file.
TimeGovernor * time_
Time is shared with flow equation.
std::vector< std::string > raw_output_strings_
Output lines of cells.
Field< 3, FieldValue< 3 >::Scalar > ref_divergence
Field< 3, FieldValue< 3 >::Scalar > cross_section
Field< 3, FieldValue< 3 >::Scalar > conductivity
Field< 3, FieldValue< 3 >::VectorFixed > field_ele_velocity
Field< 3, FieldValue< 3 >::VectorFixed > ref_velocity
Precompute l2 difference outputs.
Field< 3, FieldValue< 3 >::Scalar > ref_pressure
auto & orig_nodes_order() const
Definition: accessors.hh:240
unsigned int input_id() const
Return the element ID in the input mesh. Should be only used for special output.
Definition: accessors.hh:220
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:230
double measure() const
Computes the measure of the element.
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
Directing class of FieldValueCache.
unsigned int n_sides() const
Definition: elements.h:131
unsigned int n_nodes() const
Definition: elements.h:125
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
Definition: fe_values.hh:223
const FEValuesViews::Vector< spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
Definition: fe_values.hh:277
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
Definition: fe_values.cc:137
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:537
arma::vec::fixed< spacedim > point(const unsigned int point_no)
Return coordinates of the quadrature point in the actual cell system.
Definition: fe_values.hh:236
Raviart-Thomas element of order 0.
Definition: fe_rt.hh:61
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
@ output_file
Definition: file_path.hh:69
Generic class of assemblation.
void end() override
Implements AssemblyBase::end.
L2DifferenceAssembly(EqFields *eq_fields, EqData *eq_data)
Constructor.
arma::vec3 flux_in_q_point_
Precomputed flux in quadrature point.
FieldSet used_fields_
Sub field set contains fields used in calculation.
DarcyLMH::EqFields EqFields
std::vector< double > fluxes_
Precomputed fluxes on element sides.
DarcyFlowMHOutput::DiffEqData EqData
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void begin() override
Implements AssemblyBase::begin.
double cross_
Field results.
virtual ~L2DifferenceAssembly()
Destructor.
static constexpr const char * name()
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
FE_RT0< dim > fe_rt_
FEValues for velocity.
EqFields * eq_fields_
Data objects shared with Flow equation.
std::shared_ptr< FE_P_disc< dim > > fe_p0_
arma::vec3 q_point_
Local coords of quadrature point.
arma::vec3 ref_flux_
Field result.
FieldSet used_fields_
Sub field set contains fields used in calculation.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
DarcyFlowMHOutput::RawOutputEqData EqData
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
static constexpr const char * name()
void begin() override
Implements AssemblyBase::begin.
OutputInternalFlowAssembly(EqFields *eq_fields, EqData *eq_data)
Constructor.
virtual ~OutputInternalFlowAssembly()
Destructor.
void end() override
Implements AssemblyBase::end.
EqFields * eq_fields_
Data objects shared with Flow equation.
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:209
double t() const
Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
Output class for darcy_flow_mh model.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Definitions of Raviart-Thomas finite elements.
Class FESystem for compound finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ bulk
int IntIdx
Definition: index_types.hh:25
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
unsigned int uint
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
std::string format(CStringRef format_str, ArgList args)
Definition: format.h:3141
Definitions of particular quadrature rules on simplices.
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.
Definition: update_flags.hh:68
@ update_values
Shape function values.
Definition: update_flags.hh:87
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.