Flow123d  master-c754b67
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 
189  static const IT::Selection & get_output_selection();
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_ */
ConvectionTransport::set_balance_object
void set_balance_object(std::shared_ptr< Balance > balance) override
Definition: transport.cc:539
ConvectionTransport::is_bc_term_scaled
bool is_bc_term_scaled
Definition: transport.h:300
ConvectionTransport::EqFields::region_id
Field< 3, FieldValue< 3 >::Scalar > region_id
Definition: transport.h:108
ConvectionTransport::FactoryBaseType
ConcentrationTransportBase FactoryBaseType
Definition: transport.h:185
ConvectionTransport::EqData::bcvcorr
Vec * bcvcorr
Boundary condition correction vector.
Definition: transport.h:168
ConvectionTransport::EqData::mass_diag
Vec mass_diag
Definition: transport.h:158
ConvectionTransport::EqData::tm_diag
vector< VectorMPI > tm_diag
additions to PETSC transport matrix on the diagonal - from sources (for each substance)
Definition: transport.h:167
vector_mpi.hh
ConvectionTransport::get_input_type
static const Input::Type::Record & get_input_type()
Definition: transport.cc:73
ConvectionTransport::substances
SubstanceList & substances() override
Returns reference to the vector of substance names.
Definition: transport.h:275
Input
Abstract linear system class.
Definition: balance.hh:40
ConvectionTransport::EqData::subst_idx
vector< unsigned int > subst_idx
List of indices used to call balance methods for a set of quantities.
Definition: transport.h:156
SubstanceList::size
unsigned int size() const
Definition: substance.hh:87
ConvectionTransport::matrix_mpi_assembly_
GenericAssembly< MatrixMpiAssemblyConvection > * matrix_mpi_assembly_
Definition: transport.h:325
SubstanceList
Definition: substance.hh:70
ConcSourcesBdrAssemblyConvection
Definition: assembly_convection.hh:183
FEValues::normal_vector
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
Balance
Definition: balance.hh:119
TransportEqFields::flow_flux
Field< 3, FieldValue< 3 >::VectorFixed > flow_flux
Flow flux, can be result of water flow model.
Definition: transport_operator_splitting.hh:153
ConvectionTransport::update_after_reactions
void update_after_reactions(bool) override
Not used in this class.
Definition: transport.h:224
FEValues::JxW
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
ConvectionTransport::EqData::cfl_source_
VectorMPI cfl_source_
Parallel vector for source term contribution to CFL condition.
Definition: transport.h:164
ConvectionTransport::EqFields::output_fields
EquationOutput output_fields
Fields indended for output, i.e. all input fields plus those representing solution.
Definition: transport.h:115
ConvectionTransport::EqFields::conc_mobile
MultiField< 3, FieldValue< 3 >::Scalar > conc_mobile
Calculated concentrations in the mobile zone.
Definition: transport.h:111
EquationBase::time
TimeGovernor & time()
Definition: equation.hh:151
fe_values.hh
Class FEValues calculates finite element data on the actual cells such as shape function values,...
ConvectionTransport::EqData::is_mass_diag_changed
bool is_mass_diag_changed
Flag indicates that porosity or cross_section changed during last time.
Definition: transport.h:150
VectorMPI::resize
void resize(unsigned int local_size)
Definition: vector_mpi.cc:50
substance.hh
Classes for storing substance data.
ConvectionTransport::EqData::corr_vec
vector< VectorMPI > corr_vec
Definition: transport.h:163
ConvectionTransport::~ConvectionTransport
virtual ~ConvectionTransport()
Definition: transport.cc:206
std::vector< unsigned int >
ConvectionTransport::eq_data_
std::shared_ptr< EqData > eq_data_
Definition: transport.h:292
ConvectionTransport::EqFields::side_flux
double side_flux(SidePoint &side_p, FEValues< 3 > &fe_side_values)
Calculate flux on given side point.
Definition: transport.h:95
ConvectionTransport::n_substances
unsigned int n_substances() override
Returns number of transported substances.
Definition: transport.h:271
ConvectionTransport::EqFields::bc_conc
BCMultiField< 3, FieldValue< 3 >::Scalar > bc_conc
Definition: transport.h:103
SidePoint
General point a+ side_begin_ + ccessor allow iterate over quadrature points of given side defined in ...
Definition: eval_subset.hh:116
ConvectionTransport::vcumulative_corr
Vec * vcumulative_corr
Definition: transport.h:313
MassAssemblyConvection
Definition: assembly_convection.hh:35
ConvectionTransport::output_stream_
std::shared_ptr< OutputTime > output_stream_
Definition: transport.h:318
type_base.hh
index_types.hh
ConvectionTransport::init_cond_assembly_
GenericAssembly< InitCondAssemblyConvection > * init_cond_assembly_
Definition: transport.h:323
ConvectionTransport::EqData::alloc_transport_structs_mpi
void alloc_transport_structs_mpi(unsigned int lsize)
Definition: transport.h:134
ConvectionTransport::EqFields::init_conc
MultiField< 3, FieldValue< 3 >::Scalar > init_conc
Initial concentrations.
Definition: transport.h:106
ConvectionTransport::EqData::balance_
std::shared_ptr< Balance > balance_
object for calculation and writing the mass balance to file, shared with EquationBase.
Definition: transport.h:147
FEValues< 3 >
InitCondAssemblyConvection
Definition: assembly_convection.hh:118
equation_output.hh
Distribution
Definition: distribution.hh:50
EquationBase::balance
std::shared_ptr< Balance > balance() const
Definition: equation.hh:189
ConvectionTransport::EqData
Definition: transport.h:119
ConvectionTransport::alloc_transport_structs_mpi
void alloc_transport_structs_mpi()
Definition: transport.cc:251
ConvectionTransport::EqData::EqData
EqData()
Definition: transport.h:122
ConvectionTransport::EqData::transport_matrix_time
double transport_matrix_time
Definition: transport.h:175
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
ConvectionTransport::input_rec
const Input::Record input_rec
Record with input specification.
Definition: transport.h:316
ConvectionTransport::evaluate_time_constraint
bool evaluate_time_constraint(double &time_constraint) override
Definition: transport.cc:326
ConvectionTransport::registrar
static const int registrar
Registrar of class to factory.
Definition: transport.h:286
ConvectionTransport::EqData::max_edg_sides
unsigned int max_edg_sides
Maximal number of edge sides (evaluate from dim 1,2,3)
Definition: transport.h:180
ConvectionTransport::zero_time_step
void zero_time_step() override
Definition: transport.cc:294
quadrature_lib.hh
Definitions of particular quadrature rules on simplices.
type_generic.hh
MatrixMpiAssemblyConvection
Definition: assembly_convection.hh:358
ConvectionTransport::EqFields
Definition: transport.h:88
ConvectionTransport::EqData::n_substances
unsigned int n_substances()
Returns number of transported substances.
Definition: transport.h:126
accessors.hh
EquationOutput
Definition: equation_output.hh:46
ConvectionTransport::vpmass_diag
Vec vpmass_diag
Definition: transport.h:309
TimeGovernor
Basic time management functionality for unsteady (and steady) solvers (class Equation).
Definition: time_governor.hh:310
ConvectionTransport::output_data
virtual void output_data() override
Write computed fields.
Definition: transport.cc:527
bc_multi_field.hh
ConvectionTransport::conc_sources_bdr_assembly_
GenericAssembly< ConcSourcesBdrAssemblyConvection > * conc_sources_bdr_assembly_
Definition: transport.h:324
field_values.hh
TimeMark::Type
Definition: time_marks.hh:60
ConvectionTransport::EqFields::EqFields
EqFields()
Definition: transport.cc:90
OutputTime
The class for outputting data during time.
Definition: output_time.hh:51
ConvectionTransport
Definition: transport.h:86
ConvectionTransport::EqData::transport_bc_time
double transport_bc_time
Time of the last update of the boundary condition terms.
Definition: transport.h:169
ConvectionTransport::eq_fields_
std::shared_ptr< EqFields > eq_fields_
Definition: transport.h:291
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:65
ConvectionTransport::EqFields::conc_mobile_fe
FieldFEScalarVec conc_mobile_fe
Underlaying FieldFE for each substance of conc_mobile.
Definition: transport.h:112
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:182
ConvectionTransport::EqData::set_time_governor
void set_time_governor(TimeGovernor *time)
Definition: transport.h:129
ConvectionTransport::get_component_vec
Vec get_component_vec(unsigned int sbi) override
Return PETSc vector with solution for sbi-th substance.
Definition: transport.cc:199
ConvectionTransport::get_p0_interpolation
FieldFEScalarVec & get_p0_interpolation() override
Getter for P0 interpolation by FieldFE.
Definition: transport.cc:522
BCMultiField
Definition: bc_multi_field.hh:29
ConcentrationTransportBase
Definition: transport_operator_splitting.hh:61
ConvectionTransport::is_src_term_scaled
bool is_src_term_scaled
Definition: transport.h:300
Mesh
Definition: mesh.h:362
TransportOperatorSplitting
Coupling of a transport model with a reaction model by operator splitting.
Definition: transport_operator_splitting.hh:205
ConvectionTransport::EqFields::subdomain
Field< 3, FieldValue< 3 >::Scalar > subdomain
Definition: transport.h:109
multi_field.hh
ConvectionTransport::EqData::sources_changed_
bool sources_changed_
Flag indicates that sources part of equation was changed during last time.
Definition: transport.h:161
ConvectionTransport::target_mark_type
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
ConvectionTransport::EqData::tm
Mat tm
PETSc transport matrix.
Definition: transport.h:166
MultiField
Class for representation of a vector of fields of the same physical quantity.
Definition: multi_field.hh:87
ConvectionTransport::update_solution
void update_solution() override
Definition: transport.cc:395
region.hh
ConvectionTransport::EqData::cfl_flow_
VectorMPI cfl_flow_
Parallel vector for flow contribution to CFL condition.
Definition: transport.h:165
ConvectionTransport::EqFields::~EqFields
virtual ~EqFields()
Definition: transport.h:92
ConvectionTransport::get_subst_idx
const vector< unsigned int > & get_subst_idx() override
Return substance indices used in balance.
Definition: transport.h:249
ConvectionTransport::EqData::substances_
SubstanceList substances_
Transported substances.
Definition: transport.h:153
ConvectionTransport::ConvectionTransport
ConvectionTransport(Mesh &init_mesh, const Input::Record in_rec)
Definition: transport.cc:124
ConvectionTransport::mass_assembly_
GenericAssembly< MassAssemblyConvection > * mass_assembly_
general assembly objects, hold assembly objects of appropriate dimension
Definition: transport.h:322
ConvectionTransport::set_target_time
void set_target_time(double target_time) override
Definition: transport.cc:499
ConvectionTransport::get_output_selection
static const IT::Selection & get_output_selection()
ConvectionTransport::set_output_stream
void set_output_stream(std::shared_ptr< OutputTime > stream) override
Setter for output stream.
Definition: transport.h:257
time_marks.hh
ConvectionTransport::alloc_transport_vectors
void alloc_transport_vectors()
Definition: transport.cc:243
ConvectionTransport::EqData::time_
TimeGovernor * time_
Definition: transport.h:170
ConvectionTransport::EqData::~EqData
virtual ~EqData()
Definition: transport.h:123
ConvectionTransport::EqData::is_convection_matrix_scaled
bool is_convection_matrix_scaled
Flag indicates the state of object.
Definition: transport.h:177
VectorMPI
Definition: vector_mpi.hh:43
ConvectionTransport::vconc_out_scatter
VecScatter vconc_out_scatter
Definition: transport.h:308
Field
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
TransportEqFields
Definition: transport_operator_splitting.hh:137
ASSERT_PTR
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
ConvectionTransport::vpconc
Vec * vpconc
Definition: transport.h:312
ConvectionTransport::cfl_max_step
double cfl_max_step
Time step constraint coming from CFL condition.
Definition: transport.h:305
ConvectionTransport::initialize
void initialize() override
Definition: transport.cc:148
ConvectionTransport::EqData::dh_
std::shared_ptr< DOFHandlerMultiDim > dh_
Definition: transport.h:144
GenericAssembly
Generic class of assemblation.
Definition: generic_assembly.hh:165
transport_operator_splitting.hh
field.hh
ConcentrationTransportBase::FieldFEScalarVec
std::vector< std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > > FieldFEScalarVec
Definition: transport_operator_splitting.hh:64
ConvectionTransport::compute_p0_interpolation
void compute_p0_interpolation() override
Definition: transport.h:221
IntDim
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19