Flow123d  master-f44eb46
transport.h
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 transport.h
15  * @brief
16  * @todo
17  * - remove transport_sources
18  * - in create_transport_matric_mpi, there there is condition edge_flow > ZERO
19  * this makes matrix sparser, but can lead to elements without outflow and other problems
20  * when there are big differences in fluxes, more over it doesn't work if overall flow is very small
21  */
22 
23 #ifndef TRANSPORT_H_
24 #define TRANSPORT_H_
25 
26 
27 #include <memory> // for shared_ptr
28 #include <vector> // for vector
29 #include <petscmat.h>
30 #include "fem/fe_values.hh" // for FEValues
31 #include "fields/field.hh" // for Field
32 #include "fields/bc_multi_field.hh"
33 #include "fields/field_values.hh"
34 #include "fields/multi_field.hh"
35 #include "la/vector_mpi.hh"
37 #include "input/type_base.hh" // for Array
38 #include "input/type_generic.hh" // for Instance
39 #include "input/accessors.hh"
40 #include "system/index_types.hh"
41 #include "mesh/region.hh" // for RegionSet
42 #include "petscvec.h" // for Vec, _p_Vec
43 #include "tools/time_marks.hh" // for TimeMark, TimeM...
44 #include "transport/substance.hh" // for SubstanceList
47 #include "la/vector_mpi.hh"
48 
49 class OutputTime;
50 class Mesh;
51 class Distribution;
52 class Balance;
53 namespace Input {
54  namespace Type {
55  class Record;
56  class Selection;
57  }
58 }
59 template<unsigned int dim> class MassAssemblyConvection;
60 template<unsigned int dim> class InitCondAssemblyConvection;
61 template<unsigned int dim> class ConcSourcesBdrAssemblyConvection;
62 template<unsigned int dim> class MatrixMpiAssemblyConvection;
63 template< template<IntDim...> class DimAssembly> class GenericAssembly;
64 
65 
66 //=============================================================================
67 // TRANSPORT
68 //=============================================================================
69 
70 /**
71  * TODO:
72  * - doxy documentation
73  * - make separate method for changing time step (rescaling only), reassembly matrix only when data are changed
74  *
75  * - needs support from EqData to determine next change of the data for 1) transport matrix 2) source vector
76  * this allows us precisely choose interval where to fix timestep
77  * : field->next_change_time() - returns time jump or actual time in case of time dep. field
78  */
79 
80 
81 
82 /**
83  * Class that implements explicit finite volumes scheme with upwind. The timestep is given by CFL condition.
84  *
85  */
87 public:
88  class EqFields : public TransportEqFields {
89  public:
90 
91  EqFields();
92  virtual ~EqFields() {};
93 
94  /// Calculate flux on given side point.
95  inline double side_flux(SidePoint &side_p, FEValues<3> &fe_side_values) {
96  return arma::dot(this->flow_flux(side_p), fe_side_values.normal_vector(0)) * fe_side_values.JxW(0);
97  }
98 
99  /**
100  * Boundary conditions (Dirichlet) for concentrations.
101  * They are applied only on water inflow part of the boundary.
102  */
104 
105  /// Initial concentrations.
107 
110 
111  MultiField<3, FieldValue<3>::Scalar> conc_mobile; ///< Calculated concentrations in the mobile zone.
112  FieldFEScalarVec conc_mobile_fe; ///< Underlaying FieldFE for each substance of conc_mobile.
113 
114  /// Fields indended for output, i.e. all input fields plus those representing solution.
116  };
117 
118 
119  class EqData {
120  public:
121 
122  EqData() : is_mass_diag_changed(false), cfl_source_(PETSC_COMM_WORLD), cfl_flow_(PETSC_COMM_WORLD) {}
123  virtual ~EqData() {};
124 
125  /// Returns number of transported substances.
126  inline unsigned int n_substances()
127  { return substances_.size(); }
128 
130  ASSERT_PTR(time);
131  this->time_ = time;
132  }
133 
134  void alloc_transport_structs_mpi(unsigned int lsize) {
135  this->cfl_flow_.resize(lsize);
136  this->cfl_source_.resize(lsize);
137  }
138 
139  /**
140  * Temporary solution how to pass velocity field form the flow model.
141  * TODO: introduce FieldDiscrete -containing true DOFHandler and data vector and pass such object together with other
142  * data. Possibly make more general set_data method, allowing setting data given by name. needs support from EqDataBase.
143  */
144  std::shared_ptr<DOFHandlerMultiDim> dh_;
145 
146  /// object for calculation and writing the mass balance to file, shared with EquationBase.
147  std::shared_ptr<Balance> balance_;
148 
149  /// Flag indicates that porosity or cross_section changed during last time.
151 
152  /// Transported substances.
154 
155  /// List of indices used to call balance methods for a set of quantities.
157 
158  Vec mass_diag; // diagonal entries in pass matrix (cross_section * porosity)
159 
160  /// Flag indicates that sources part of equation was changed during last time.
162 
164  VectorMPI cfl_source_; ///< Parallel vector for source term contribution to CFL condition.
165  VectorMPI cfl_flow_; ///< Parallel vector for flow contribution to CFL condition.
166  Mat tm; ///< PETSc transport matrix
167  vector<VectorMPI> tm_diag; ///< additions to PETSC transport matrix on the diagonal - from sources (for each substance)
168  Vec *bcvcorr; ///< Boundary condition correction vector
169  double transport_bc_time; ///< Time of the last update of the boundary condition terms.
171 
172  /// Time when the transport matrix was created.
173  /// TODO: when we have our own classes for LA objects, we can use lazy dependence to check
174  /// necessity for matrix update
176 
177  bool is_convection_matrix_scaled; ///< Flag indicates the state of object
178 
179  /// Maximal number of edge sides (evaluate from dim 1,2,3)
180  unsigned int max_edg_sides;
181 
182  };
183 
184 
186 
187  static const Input::Type::Record & get_input_type();
188 
190 
191  /**
192  * Constructor.
193  */
194  ConvectionTransport(Mesh &init_mesh, const Input::Record in_rec);
195  /**
196  * TODO: destructor
197  */
198  virtual ~ConvectionTransport();
199 
200  void initialize() override;
201 
202  /**
203  * Initialize solution at zero time.
204  */
205  void zero_time_step() override;
206  /**
207  * Evaluates CFL condition.
208  * Assembles the transport matrix and vector (including sources, bc terms).
209  * @param time_constraint is the value CFL constraint (return parameter)
210  * @return true if CFL is changed since previous step, false otherwise
211  */
212  bool evaluate_time_constraint(double &time_constraint) override;
213  /**
214  * Calculates one time step of explicit transport.
215  */
216  void update_solution() override;
217 
218  /** Compute P0 interpolation of the solution (used reaction term).
219  * Empty - solution is already P0 interpolation.
220  */
221  void compute_p0_interpolation() override {};
222 
223  /// Not used in this class.
224  void update_after_reactions(bool) override {};
225 
226  /**
227  * Set time interval which is considered as one time step by TransportOperatorSplitting.
228  * In particular the velocity field doesn't change over this interval.
229  *
230  * Dependencies:
231  *
232  * velocity, porosity -> matrix, source_vector
233  * matrix -> time_step
234  *
235  * data_read_times -> time_step (not necessary if we won't stick to jump times)
236  * data -> source_vector
237  * time_step -> scaling
238  *
239  *
240  *
241  */
242  void set_target_time(double target_time) override;
243 
244  /**
245  * Use Balance object from upstream equation (e.g. in various couplings) instead of own instance.
246  */
247  void set_balance_object(std::shared_ptr<Balance> balance) override;
248 
250  { return eq_data_->subst_idx; }
251 
252  /**
253  * @brief Write computed fields.
254  */
255  virtual void output_data() override;
256 
257  void set_output_stream(std::shared_ptr<OutputTime> stream) override
258  { output_stream_ = stream; }
259 
260 
261  /**
262  * Getters.
263  */
264 
265  /// Getter for P0 interpolation by FieldFE.
267 
268  Vec get_component_vec(unsigned int sbi) override;
269 
270  /// Returns number of transported substances.
271  inline unsigned int n_substances() override
272  { return eq_data_->substances_.size(); }
273 
274  /// Returns reference to the vector of substance names.
275  inline SubstanceList &substances() override
276  { return eq_data_->substances_; }
277 
278 private:
279 
282 
283 
284 
285  /// Registrar of class to factory
286  static const int registrar;
287 
288  /**
289  * Parameters of the equation, some are shared with other implementations since EqFields is derived from TransportBase::TransportEqFields
290  */
291  std::shared_ptr<EqFields> eq_fields_;
292  std::shared_ptr<EqData> eq_data_;
293 
294  //@{
295  /**
296  * Flag indicates the state of object (transport matrix or source or boundary term).
297  * If false, the object is freshly assembled and not rescaled.
298  * If true, the object is scaled (not necessarily with the current time step).
299  */
301 
302  //@}
303 
304  TimeMark::Type target_mark_type; ///< TimeMark type for time marks denoting end of every time interval where transport matrix remains constant.
305  double cfl_max_step; ///< Time step constraint coming from CFL condition.
306 
307 
308  VecScatter vconc_out_scatter;
309  Vec vpmass_diag; // diagonal entries in mass matrix from last time (cross_section * porosity)
310 
311  ///
312  Vec *vpconc; // previous concentration vector
314 
315  /// Record with input specification.
317 
318  std::shared_ptr<OutputTime> output_stream_;
319 
320 
321  /// general assembly objects, hold assembly objects of appropriate dimension
326 
328 };
329 #endif /* TRANSPORT_H_ */
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
std::vector< std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > > FieldFEScalarVec
bool sources_changed_
Flag indicates that sources part of equation was changed during last time.
Definition: transport.h:161
bool is_convection_matrix_scaled
Flag indicates the state of object.
Definition: transport.h:177
void alloc_transport_structs_mpi(unsigned int lsize)
Definition: transport.h:134
Vec * bcvcorr
Boundary condition correction vector.
Definition: transport.h:168
VectorMPI cfl_source_
Parallel vector for source term contribution to CFL condition.
Definition: transport.h:164
std::shared_ptr< Balance > balance_
object for calculation and writing the mass balance to file, shared with EquationBase.
Definition: transport.h:147
VectorMPI cfl_flow_
Parallel vector for flow contribution to CFL condition.
Definition: transport.h:165
unsigned int max_edg_sides
Maximal number of edge sides (evaluate from dim 1,2,3)
Definition: transport.h:180
double transport_bc_time
Time of the last update of the boundary condition terms.
Definition: transport.h:169
void set_time_governor(TimeGovernor *time)
Definition: transport.h:129
bool is_mass_diag_changed
Flag indicates that porosity or cross_section changed during last time.
Definition: transport.h:150
unsigned int n_substances()
Returns number of transported substances.
Definition: transport.h:126
vector< VectorMPI > tm_diag
additions to PETSC transport matrix on the diagonal - from sources (for each substance)
Definition: transport.h:167
Mat tm
PETSc transport matrix.
Definition: transport.h:166
vector< VectorMPI > corr_vec
Definition: transport.h:163
SubstanceList substances_
Transported substances.
Definition: transport.h:153
vector< unsigned int > subst_idx
List of indices used to call balance methods for a set of quantities.
Definition: transport.h:156
std::shared_ptr< DOFHandlerMultiDim > dh_
Definition: transport.h:144
Field< 3, FieldValue< 3 >::Scalar > subdomain
Definition: transport.h:109
MultiField< 3, FieldValue< 3 >::Scalar > init_conc
Initial concentrations.
Definition: transport.h:106
Field< 3, FieldValue< 3 >::Scalar > region_id
Definition: transport.h:108
FieldFEScalarVec conc_mobile_fe
Underlaying FieldFE for each substance of conc_mobile.
Definition: transport.h:112
double side_flux(SidePoint &side_p, FEValues< 3 > &fe_side_values)
Calculate flux on given side point.
Definition: transport.h:95
EquationOutput output_fields
Fields indended for output, i.e. all input fields plus those representing solution.
Definition: transport.h:115
BCMultiField< 3, FieldValue< 3 >::Scalar > bc_conc
Definition: transport.h:103
MultiField< 3, FieldValue< 3 >::Scalar > conc_mobile
Calculated concentrations in the mobile zone.
Definition: transport.h:111
ConcentrationTransportBase FactoryBaseType
Definition: transport.h:185
void set_output_stream(std::shared_ptr< OutputTime > stream) override
Setter for output stream.
Definition: transport.h:257
SubstanceList & substances() override
Returns reference to the vector of substance names.
Definition: transport.h:275
static const IT::Selection & get_output_selection()
FieldFEScalarVec & get_p0_interpolation() override
Getter for P0 interpolation by FieldFE.
Definition: transport.cc:522
unsigned int n_substances() override
Returns number of transported substances.
Definition: transport.h:271
VecScatter vconc_out_scatter
Definition: transport.h:308
void initialize() override
Definition: transport.cc:148
void set_target_time(double target_time) override
Definition: transport.cc:499
void update_solution() override
Definition: transport.cc:395
GenericAssembly< InitCondAssemblyConvection > * init_cond_assembly_
Definition: transport.h:323
void alloc_transport_structs_mpi()
Definition: transport.cc:251
std::shared_ptr< EqFields > eq_fields_
Definition: transport.h:291
void zero_time_step() override
Definition: transport.cc:294
GenericAssembly< MassAssemblyConvection > * mass_assembly_
general assembly objects, hold assembly objects of appropriate dimension
Definition: transport.h:322
void alloc_transport_vectors()
Definition: transport.cc:243
GenericAssembly< MatrixMpiAssemblyConvection > * matrix_mpi_assembly_
Definition: transport.h:325
void update_after_reactions(bool) override
Not used in this class.
Definition: transport.h:224
ConvectionTransport(Mesh &init_mesh, const Input::Record in_rec)
Definition: transport.cc:124
const Input::Record input_rec
Record with input specification.
Definition: transport.h:316
virtual void output_data() override
Write computed fields.
Definition: transport.cc:527
Vec get_component_vec(unsigned int sbi) override
Return PETSc vector with solution for sbi-th substance.
Definition: transport.cc:199
std::shared_ptr< OutputTime > output_stream_
Definition: transport.h:318
void set_balance_object(std::shared_ptr< Balance > balance) override
Definition: transport.cc:539
static const Input::Type::Record & get_input_type()
Definition: transport.cc:73
TimeMark::Type target_mark_type
TimeMark type for time marks denoting end of every time interval where transport matrix remains const...
Definition: transport.h:304
void compute_p0_interpolation() override
Definition: transport.h:221
GenericAssembly< ConcSourcesBdrAssemblyConvection > * conc_sources_bdr_assembly_
Definition: transport.h:324
bool evaluate_time_constraint(double &time_constraint) override
Definition: transport.cc:326
double cfl_max_step
Time step constraint coming from CFL condition.
Definition: transport.h:305
static const int registrar
Registrar of class to factory.
Definition: transport.h:286
std::shared_ptr< EqData > eq_data_
Definition: transport.h:292
const vector< unsigned int > & get_subst_idx() override
Return substance indices used in balance.
Definition: transport.h:249
virtual ~ConvectionTransport()
Definition: transport.cc:206
std::shared_ptr< Balance > balance() const
Definition: equation.hh:189
TimeGovernor & time()
Definition: equation.hh:151
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
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
Definition: fe_values.hh:257
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
Generic class of assemblation.
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Record type proxy class.
Definition: type_record.hh:182
Template for classes storing finite set of named values.
Definition: mesh.h:362
Class for representation of a vector of fields of the same physical quantity.
Definition: multi_field.hh:87
The class for outputting data during time.
Definition: output_time.hh:51
General point a+ side_begin_ + ccessor allow iterate over quadrature points of given side defined in ...
Definition: eval_subset.hh:116
unsigned int size() const
Definition: substance.hh:87
Basic time management functionality for unsteady (and steady) solvers (class Equation).
Field< 3, FieldValue< 3 >::VectorFixed > flow_flux
Flow flux, can be result of water flow model.
Coupling of a transport model with a reaction model by operator splitting.
void resize(unsigned int local_size)
Definition: vector_mpi.cc:50
Class FEValues calculates finite element data on the actual cells such as shape function values,...
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Abstract linear system class.
Definition: balance.hh:40
Definitions of particular quadrature rules on simplices.
Classes for storing substance data.