Flow123d  DF_patch_fe_data_tables-b828b90
concentration_model.cc
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 concentration_model.cc
15  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
16  * @author Jan Stebel
17  */
18 
19 #include "input/input_type.hh"
20 #include "mesh/mesh.h"
21 #include "mesh/accessors.hh"
22 
24 #include "concentration_model.hh"
25 #include "tools/unit_si.hh"
26 #include "coupling/balance.hh"
27 #include "fields/field_model.hh"
28 
29 
30 
31 using namespace std;
32 using namespace Input::Type;
33 
34 
35 /*******************************************************************************
36  * Functors of FieldModels
37  */
38 using Sclr = double;
39 using Vect = arma::vec3;
40 using Tens = arma::mat33;
41 
42 // Functor computing velocity norm
44  inline Sclr operator() (Vect vel) {
45  return arma::norm(vel, 2);
46  }
47 };
48 
49 // Functor computing mass matrix coefficients (cross_section * water_content)
51  inline Sclr operator() (Sclr csec, Sclr wcont) {
52  return csec * wcont;
53  }
54 };
55 
56 // Functor computing retardation coefficients:
57 // (1-porosity) * rock_density * sorption_coefficient * rock_density
59  inline Sclr operator() (Sclr csec, Sclr por_m, Sclr rho_s, Sclr sorp_mult) {
60  return (1.-por_m)*rho_s*sorp_mult*csec;
61  }
62 };
63 
64 // Functor computing sources density output (cross_section * sources_density)
66  inline Sclr operator() (Sclr csec, Sclr sdens) {
67  return csec * sdens;
68  }
69 };
70 
71 // Functor computing sources sigma output (cross_section * sources_sigma)
73  inline Sclr operator() (Sclr csec, Sclr ssigma) {
74  return csec * ssigma;
75  }
76 };
77 
78 // Functor computing sources concentration output (sources_conc)
80  inline Sclr operator() (Sclr sconc) {
81  return sconc;
82  }
83 };
84 
85 // Functor computing advection coefficient (velocity)
87  inline Vect operator() (Vect velocity) {
88  return velocity;
89  }
90 };
91 
92 // Functor computing diffusion coefficient (see notes in function)
94  inline Tens operator() (Tens diff_m, Vect velocity, Sclr v_norm, Sclr alphaL, Sclr alphaT, Sclr water_content, Sclr porosity, Sclr c_sec) {
95 
96  // used tortuosity model dues to Millington and Quirk(1961) (should it be with power 10/3 ?)
97  // for an overview of other models see: Chou, Wu, Zeng, Chang (2011)
98  double tortuosity = pow(water_content, 7.0 / 3.0)/ (porosity * porosity);
99 
100  // result
101  Tens K;
102 
103  // Note that the velocity vector is in fact the Darcian flux,
104  // so we need not to multiply vnorm by water_content and cross_section.
105  //K = ((alphaL-alphaT) / vnorm) * K + (alphaT*vnorm + Dm*tortuosity*cross_cut*water_content) * arma::eye(3,3);
106 
107  if (fabs(v_norm) > 0) {
108  /*
109  for (int i=0; i<3; i++)
110  for (int j=0; j<3; j++)
111  K(i,j) = (velocity[i]/vnorm)*(velocity[j]);
112  */
113  K = ((alphaL - alphaT) / v_norm) * arma::kron(velocity.t(), velocity);
114 
115  //arma::mat33 abs_diff_mat = arma::abs(K - kk);
116  //double diff = arma::min( arma::min(abs_diff_mat) );
117  //ASSERT_PERMANENT( diff < 1e-12 )(diff)(K)(kk);
118  } else
119  K.zeros();
120 
121  // Note that the velocity vector is in fact the Darcian flux,
122  // so to obtain |v| we have to divide vnorm by porosity and cross_section.
123  K += alphaT*v_norm*arma::eye(3,3) + diff_m*(tortuosity*c_sec*water_content);
124 
125  return K;
126  }
127 };
128 
129 
130 
131 
132 
135 {
136  *this+=bc_type
137  .name("bc_type")
138  .description(
139  "Type of boundary condition.")
141  .input_default("\"inflow\"")
144  *this+=bc_dirichlet_value
145  .name("bc_conc")
146  .units( UnitSI().kg().m(-3) )
147  .description("Dirichlet boundary condition (for each substance).")
148  .input_default("0.0")
149  .flags_add( in_rhs );
150  *this+=bc_flux
152  .name("bc_flux")
153  .description("Flux in Neumann boundary condition.")
154  .units( UnitSI().kg().m().s(-1).md() )
155  .input_default("0.0")
156  .flags_add(FieldFlag::in_rhs);
157  *this+=bc_robin_sigma
159  .name("bc_robin_sigma")
160  .description("Conductivity coefficient in Robin boundary condition.")
161  .units( UnitSI().m(4).s(-1).md() )
162  .input_default("0.0")
164  *this+=init_condition
165  .name("init_conc")
166  .units( UnitSI().kg().m(-3) )
167  .description("Initial values for concentration of substances.")
168  .input_default("0.0");
169  *this+=disp_l
170  .name("disp_l")
171  .description("Longitudinal dispersivity in the liquid (for each substance).")
172  .units( UnitSI().m() )
173  .input_default("0.0")
175  *this+=disp_t
176  .name("disp_t")
177  .description("Transverse dispersivity in the liquid (for each substance).")
178  .units( UnitSI().m() )
179  .input_default("0.0")
181  *this+=diff_m
182  .name("diff_m")
183  .description("Molecular diffusivity in the liquid (for each substance).")
184  .units( UnitSI().m(2).s(-1) )
185  .input_default("0.0")
187  *this+=rock_density
188  .name("rock_density")
189  .description("Rock matrix density.")
190  .units(UnitSI().kg().m(-3))
191  .input_default("0.0")
193  *this+=sorption_coefficient
194  .name("sorption_coefficient")
195  .description("Coefficient of linear sorption.")
196  .units(UnitSI().m(3).kg(-1))
197  .input_default("0.0")
199 
200  *this+=output_field
201  .name("conc")
202  .description("Concentration solution.")
203  .units( UnitSI().kg().m(-3) )
205 
206 
207  // initiaization of FieldModels
208  *this += v_norm.name("v_norm")
209  .description("Velocity norm field.")
210  .input_default("0.0")
211  .units( UnitSI().m().s(-1) );
212 
213  *this += mass_matrix_coef.name("mass_matrix_coef")
214  .description("Matrix coefficients computed by model in mass assemblation.")
215  .input_default("0.0")
216  .units( UnitSI().m(3).md() );
217 
218  *this += retardation_coef.name("retardation_coef")
219  .description("Retardation coefficients computed by model in mass assemblation.")
220  .input_default("0.0")
221  .units( UnitSI().m(3).md() );
222  *this += sources_density_out.name("sources_density_out")
223  .description("Concentration sources output - density of substance source, only positive part is used..")
224  .input_default("0.0")
225  .units( UnitSI().kg().s(-1).md() );
226 
227  *this += sources_sigma_out.name("sources_sigma_out")
228  .description("Concentration sources - Robin type, in_flux = sources_sigma * (sources_conc - mobile_conc).")
229  .input_default("0.0")
230  .units( UnitSI().s(-1).m(3).md() );
231 
232  *this += sources_conc_out.name("sources_conc_out")
233  .description("Concentration sources output.")
234  .input_default("0.0")
235  .units( UnitSI().kg().m(-3) );
236 
237  *this += advection_coef.name("advection_coef")
238  .description("Advection coefficients model.")
239  .input_default("0.0")
240  .units( UnitSI().m().s(-1) );
241 
242  *this += diffusion_coef.name("diffusion_coef")
243  .description("Diffusion coefficients model.")
244  .input_default("0.0")
245  .units( UnitSI().m(2).s(-1) );
246 
247 }
248 
249 
251 {
252  // initialize multifield components
253  sorption_coefficient.setup_components();
254  sources_conc.setup_components();
255  sources_density.setup_components();
256  sources_sigma.setup_components();
257  diff_m.setup_components();
258  disp_l.setup_components();
259  disp_t.setup_components();
260 
261  // create FieldModels
262  v_norm.set(Model<3, FieldValue<3>::Scalar>::create(fn_conc_v_norm(), flow_flux), 0.0);
263  mass_matrix_coef.set(Model<3, FieldValue<3>::Scalar>::create(fn_conc_mass_matrix(), cross_section, water_content), 0.0);
264  retardation_coef.set(
265  Model<3, FieldValue<3>::Scalar>::create_multi(fn_conc_retardation(), cross_section, porosity, rock_density, sorption_coefficient),
266  0.0
267  );
268  sources_density_out.set(Model<3, FieldValue<3>::Scalar>::create_multi(fn_conc_sources_dens(), cross_section, sources_density), 0.0);
269  sources_sigma_out.set(Model<3, FieldValue<3>::Scalar>::create_multi(fn_conc_sources_sigma(), cross_section, sources_sigma), 0.0);
270  sources_conc_out.set(Model<3, FieldValue<3>::Scalar>::create_multi(fn_conc_sources_conc(), sources_conc), 0.0);
271  std::vector<typename Field<3, FieldValue<3>::VectorFixed>::FieldBasePtr> ad_coef_ptr_vec;
272  for (unsigned int sbi=0; sbi<sorption_coefficient.size(); sbi++)
273  ad_coef_ptr_vec.push_back( Model<3, FieldValue<3>::VectorFixed>::create(fn_conc_ad_coef(), flow_flux) );
274  advection_coef.set(ad_coef_ptr_vec, 0.0);
275  diffusion_coef.set(
276  Model<3, FieldValue<3>::TensorFixed>::create_multi(
277  fn_conc_diff_coef(), diff_m, flow_flux, v_norm, disp_l, disp_t, water_content, porosity, cross_section
278  ),
279  0.0
280  );
281 }
282 
283 
284 
286  return Selection("Solute_AdvectionDiffusion_BC_Type", "Types of boundary conditions for advection-diffusion solute transport model.")
287  .add_value(bc_inflow, "inflow",
288  "Default transport boundary condition.\n"
289  "On water inflow (($(q_w \\le 0)$)), total flux is given by the reference concentration 'bc_conc'. "
290  "On water outflow we prescribe zero diffusive flux, "
291  "i.e. the mass flows out only due to advection.")
292  .add_value(bc_dirichlet, "dirichlet",
293  "Dirichlet boundary condition (($ c = c_D $)).\n"
294  "The prescribed concentration (($c_D$)) is specified by the field 'bc_conc'.")
295  .add_value(bc_total_flux, "total_flux",
296  "Total mass flux boundary condition.\n"
297  "The prescribed total incoming flux can have the general form (($\\delta(f_N+\\sigma_R(c_R-c) )$)), "
298  "where the absolute flux (($f_N$)) is specified by the field 'bc_flux', "
299  "the transition parameter (($\\sigma_R$)) by 'bc_robin_sigma', "
300  "and the reference concentration (($c_R$)) by 'bc_conc'.")
301  .add_value(bc_diffusive_flux, "diffusive_flux",
302  "Diffusive flux boundary condition.\n"
303  "The prescribed incoming mass flux due to diffusion can have the general form (($\\delta(f_N+\\sigma_R(c_R-c) )$)), "
304  "where the absolute flux (($f_N$)) is specified by the field 'bc_flux', "
305  "the transition parameter (($\\sigma_R$)) by 'bc_robin_sigma', "
306  "and the reference concentration (($c_R$)) by 'bc_conc'.")
307  .close();
308 }
309 
310 
311 
312 
314 {
315  // Return empty selection just to provide model specific selection name and description.
316  // The fields are added by TransportDG using an auxiliary selection.
317  return IT::Selection(
318  std::string(ModelEqData::name()) + "_DG_output_fields",
319  "Selection of output fields for Diffusive Solute Transport DG model.");
320 }
321 
322 
323 IT::Record ConcentrationTransportModel::get_input_type(const string &implementation, const string &description)
324 {
325  return IT::Record(
326  std::string(ModelEqData::name()) + "_" + implementation,
327  description + " for solute transport.")
329  .declare_key("solvent_density", IT::Double(0), IT::Default("1.0"),
330  "Density of the solvent [ (($kg.m^{-3}$)) ].");
331 }
332 
333 
336 {}
337 
338 
340 {
341  solvent_density_ = in_rec.val<double>("solvent_density");
342 }
343 
344 
346 {}
347 
348 
349 void ConcentrationTransportModel::set_balance_object(std::shared_ptr<Balance> balance)
350 {
351  balance_ = balance;
352  eq_data().subst_idx_ = balance_->add_quantities(eq_data().substances_.names());
353 }
354 
355 
357 
358 
359 
static Input::Type::Abstract & get_input_type()
Common specification of the input record for secondary equations.
static constexpr const char * name()
vector< unsigned int > subst_idx_
List of indices used to call balance methods for a set of quantities.
Field< 3, FieldValue< 3 >::Scalar > rock_density
Rock matrix density.
MultiField< 3, FieldValue< 3 >::Scalar > sources_conc_out
Concentration sources - concentration output.
BCMultiField< 3, FieldValue< 3 >::Enum > bc_type
Type of boundary condition (see also BC_Type)
MultiField< 3, FieldValue< 3 >::Scalar > init_condition
Initial concentrations.
Field< 3, FieldValue< 3 >::Scalar > mass_matrix_coef
Field represents coefficients of mass matrix.
MultiField< 3, FieldValue< 3 >::Scalar > sorption_coefficient
Coefficient of linear sorption.
MultiField< 3, FieldValue< 3 >::Scalar > sources_sigma_out
Concentration sources - sigma output.
BCMultiField< 3, FieldValue< 3 >::Scalar > bc_flux
Flux value in total/diffusive flux b.c.
Field< 3, FieldValue< 3 >::Scalar > v_norm
Velocity norm field.
MultiField< 3, FieldValue< 3 >::TensorFixed > diff_m
Molecular diffusivity (for each substance).
BCMultiField< 3, FieldValue< 3 >::Scalar > bc_dirichlet_value
Prescribed concentration for Dirichlet/reference concentration for flux b.c.
MultiField< 3, FieldValue< 3 >::Scalar > output_field
MultiField< 3, FieldValue< 3 >::Scalar > retardation_coef
Field represents retardation coefficients due to sorption.
MultiField< 3, FieldValue< 3 >::Scalar > disp_t
Transversal dispersivity (for each substance).
MultiField< 3, FieldValue< 3 >::TensorFixed > diffusion_coef
Diffusion coefficients.
static const Input::Type::Selection & get_bc_type_selection()
MultiField< 3, FieldValue< 3 >::VectorFixed > advection_coef
Advection coefficients.
MultiField< 3, FieldValue< 3 >::Scalar > sources_density_out
Concentration sources - density output.
MultiField< 3, FieldValue< 3 >::Scalar > disp_l
Longitudal dispersivity (for each substance).
BCMultiField< 3, FieldValue< 3 >::Scalar > bc_robin_sigma
Transition coefficient in total/diffusive flux b.c.
virtual ModelEqData & eq_data()=0
Derived class should implement getter for ModelEqData instance.
void init_from_input(const Input::Record &in_rec) override
Read necessary data from input record.
void set_balance_object(std::shared_ptr< Balance > balance) override
void init_balance(const Input::Record &in_rec)
double solvent_density_
Density of liquid (a global constant).
ConcentrationTransportModel(Mesh &mesh, const Input::Record &in_rec)
std::shared_ptr< Balance > balance() const
Definition: equation.hh:189
Mesh & mesh()
Definition: equation.hh:181
std::shared_ptr< Balance > balance_
object for calculation and writing the mass balance to file.
Definition: equation.hh:252
FieldCommon & input_selection(Input::Type::Selection element_selection)
FieldCommon & description(const string &description)
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
static constexpr Mask in_time_term
A field is part of time term of the equation.
Definition: field_flag.hh:47
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
Definition: field_flag.hh:49
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Record type proxy class.
Definition: type_record.hh:182
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
Definition: mesh.h:362
auto disable_where(const MultiField< spacedim, typename FieldValue< spacedim >::Enum > &control_field, const vector< FieldEnum > &value_list) -> MultiField &
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
Discontinuous Galerkin method for equation of transport with dispersion.
double Sclr
#define FMT_UNUSED
Definition: posix.h:75