Flow123d  3.9.1-8f37ba8
darcy_flow_mh.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 darcy_flow_mh.cc
15  * @ingroup flow
16  * @brief Setup and solve linear system of mixed-hybrid discretization of the linear
17  * porous media flow with possible preferential flow in fractures and chanels.
18  */
19 
20 //#include <limits>
21 #include <vector>
22 //#include <iostream>
23 //#include <iterator>
24 //#include <algorithm>
25 #include <armadillo>
26 
27 #include "petscmat.h"
28 #include "petscviewer.h"
29 #include "petscerror.h"
30 #include "mpi.h"
31 
32 #include "system/system.hh"
33 #include "system/sys_profiler.hh"
34 #include "system/index_types.hh"
35 #include "input/factory.hh"
36 
37 #include "mesh/mesh.h"
38 #include "mesh/bc_mesh.hh"
39 #include "mesh/partitioning.hh"
40 #include "mesh/accessors.hh"
41 #include "mesh/range_wrapper.hh"
42 #include "la/distribution.hh"
43 #include "la/linsys.hh"
44 #include "la/linsys_PETSC.hh"
45 #include "la/linsys_BDDC.hh"
46 #include "la/schur.hh"
47 //#include "la/sparse_graph.hh"
49 #include "la/vector_mpi.hh"
50 
51 #include "flow/assembly_mh_old.hh"
52 #include "flow/darcy_flow_mh.hh"
54 #include "flow/assembly_models.hh"
55 
56 #include "tools/time_governor.hh"
58 #include "fields/field.hh"
59 #include "fields/field_values.hh"
61 #include "fields/field_fe.hh"
62 #include "fields/field_model.hh"
63 #include "fields/field_constant.hh"
64 
65 #include "coupling/balance.hh"
66 
69 
70 #include "fem/fe_p.hh"
71 
72 
73 FLOW123D_FORCE_LINK_IN_CHILD(darcy_flow_mh)
74 
75 
76 
77 
78 
79 namespace it = Input::Type;
80 
82  return it::Selection("MH_MortarMethod")
83  .add_value(NoMortar, "None", "No Mortar method is applied.")
84  .add_value(MortarP0, "P0", "Mortar space: P0 on elements of lower dimension.")
85  .add_value(MortarP1, "P1", "Mortar space: P1 on intersections, using non-conforming pressures.")
86  .close();
87 }
88 
89 
91  return it::Selection("Flow_Darcy_BC_Type")
92  .add_value(none, "none",
93  "Homogeneous Neumann boundary condition\n(zero normal flux over the boundary).")
94  .add_value(dirichlet, "dirichlet",
95  "Dirichlet boundary condition. "
96  "Specify the pressure head through the ``bc_pressure`` field "
97  "or the piezometric head through the ``bc_piezo_head`` field.")
98  .add_value(total_flux, "total_flux", "Flux boundary condition (combines Neumann and Robin type). "
99  "Water inflow equal to (($ \\delta_d(q_d^N + \\sigma_d (h_d^R - h_d) )$)). "
100  "Specify the water inflow by the ``bc_flux`` field, the transition coefficient by ``bc_robin_sigma`` "
101  "and the reference pressure head or piezometric head through ``bc_pressure`` or ``bc_piezo_head`` respectively.")
102  .add_value(seepage, "seepage",
103  "Seepage face boundary condition. Pressure and inflow bounded from above. Boundary with potential seepage flow "
104  "is described by the pair of inequalities: "
105  "(($h_d \\le h_d^D$)) and (($ -\\boldsymbol q_d\\cdot\\boldsymbol n \\le \\delta q_d^N$)), where the equality holds in at least one of them. "
106  "Caution: setting (($q_d^N$)) strictly negative "
107  "may lead to an ill posed problem since a positive outflow is enforced. "
108  "Parameters (($h_d^D$)) and (($q_d^N$)) are given by the fields ``bc_switch_pressure`` (or ``bc_switch_piezo_head``) and ``bc_flux`` respectively."
109  )
110  .add_value(river, "river",
111  "River boundary condition. For the water level above the bedrock, (($H_d > H_d^S$)), the Robin boundary condition is used with the inflow given by: "
112  "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d) )$)). For the water level under the bedrock, constant infiltration is used: "
113  "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d^S) )$)). Parameters: ``bc_pressure``, ``bc_switch_pressure``, "
114  " ``bc_sigma``, ``bc_flux``."
115  )
116  .close();
117 }
118 
120 
121  const it::Record &field_descriptor =
122  it::Record("Flow_Darcy_MH_Data",FieldCommon::field_descriptor_record_description("Flow_Darcy_MH_Data") )
123  .copy_keys( DarcyMH::EqFields().make_field_descriptor_type("Flow_Darcy_MH_Data_aux") )
124  .declare_key("bc_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
125  "Boundary piezometric head for BC types: dirichlet, robin, and river." )
126  .declare_key("bc_switch_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
127  "Boundary switch piezometric head for BC types: seepage, river." )
128  .declare_key("init_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
129  "Initial condition for the pressure given as the piezometric head." )
130  .close();
131  return field_descriptor;
132 }
133 
135 
136  it::Record ns_rec = Input::Type::Record("NonlinearSolver", "Non-linear solver settings.")
137  .declare_key("linear_solver", LinSys::get_input_type(), it::Default("{}"),
138  "Linear solver for MH problem.")
139  .declare_key("tolerance", it::Double(0.0), it::Default("1E-6"),
140  "Residual tolerance.")
141  .declare_key("min_it", it::Integer(0), it::Default("1"),
142  "Minimum number of iterations (linear solutions) to use.\nThis is usefull if the convergence criteria "
143  "does not characterize your goal well enough so it converges prematurely, possibly even without a single linear solution."
144  "If greater then 'max_it' the value is set to 'max_it'.")
145  .declare_key("max_it", it::Integer(0), it::Default("100"),
146  "Maximum number of iterations (linear solutions) of the non-linear solver.")
147  .declare_key("converge_on_stagnation", it::Bool(), it::Default("false"),
148  "If a stagnation of the nonlinear solver is detected the solver stops. "
149  "A divergence is reported by default, forcing the end of the simulation. By setting this flag to 'true', the solver "
150  "ends with convergence success on stagnation, but it reports warning about it.")
151  .close();
152 
154 
155  return it::Record("Flow_Darcy_MH", "Mixed-Hybrid solver for saturated Darcy flow.")
158  .declare_key("gravity", it::Array(it::Double(), 3,3), it::Default("[ 0, 0, -1]"),
159  "Vector of the gravity force. Dimensionless.")
161  "Input data for Darcy flow model.")
162  .declare_key("nonlinear_solver", ns_rec, it::Default("{}"),
163  "Non-linear solver for MH problem.")
164  .declare_key("output_stream", OutputTime::get_input_type(), it::Default("{}"),
165  "Output stream settings.\n Specify file format, precision etc.")
166 
167  .declare_key("output", DarcyFlowMHOutput::get_input_type(eq_fields, "Flow_Darcy_MH"),
168  IT::Default("{ \"fields\": [ \"pressure_p0\", \"velocity_p0\" ] }"),
169  "Specification of output fields and output times.")
171  "Output settings specific to Darcy flow model.\n"
172  "Includes raw output and some experimental functionality.")
173  .declare_key("balance", Balance::get_input_type(), it::Default("{}"),
174  "Settings for computing mass balance.")
175  .declare_key("n_schurs", it::Integer(0,2), it::Default("2"),
176  "Number of Schur complements to perform when solving MH system.")
177  .declare_key("mortar_method", get_mh_mortar_selection(), it::Default("\"None\""),
178  "Method for coupling Darcy flow between dimensions on incompatible meshes. [Experimental]" )
179  .close();
180 }
181 
182 
183 const int DarcyMH::registrar =
184  Input::register_class< DarcyMH, Mesh &, const Input::Record >("Flow_Darcy_MH") +
186 
187 
188 
190 {
191  *this += field_ele_pressure.name("pressure_p0")
192  .units(UnitSI().m())
194  .description("Pressure solution - P0 interpolation.");
195 
196  *this += field_edge_pressure.name("pressure_edge")
197  .units(UnitSI().m())
198  .flags(FieldFlag::input_copy)
199  .description("Pressure solution - Crouzeix-Raviart interpolation.");
200 
201  *this += field_ele_piezo_head.name("piezo_head_p0")
202  .units(UnitSI().m())
204  .description("Piezo head solution - P0 interpolation.");
205 
206  *this += field_ele_velocity.name("velocity_p0")
207  .units(UnitSI().m().s(-1))
209  .description("Velocity solution - P0 interpolation.");
210 
211  *this += flux.name("flux")
212  .units(UnitSI().m().s(-1))
214  .description("Darcy flow flux.");
215 
216  *this += anisotropy.name("anisotropy")
217  .description("Anisotropy of the conductivity tensor.")
218  .input_default("1.0")
219  .units( UnitSI::dimensionless() );
220 
221  *this += cross_section.name("cross_section")
222  .description("Complement dimension parameter (cross section for 1D, thickness for 2D).")
223  .input_default("1.0")
224  .units( UnitSI().m(3).md() );
225 
226  *this += conductivity.name("conductivity")
227  .description("Isotropic conductivity scalar.")
228  .input_default("1.0")
229  .units( UnitSI().m().s(-1) )
230  .set_limits(0.0);
231 
232  *this += sigma.name("sigma")
233  .description("Transition coefficient between dimensions.")
234  .input_default("1.0")
235  .units( UnitSI::dimensionless() );
236 
237  *this += water_source_density.name("water_source_density")
238  .description("Water source density.")
239  .input_default("0.0")
240  .units( UnitSI().s(-1) );
241 
242  *this += bc_type.name("bc_type")
243  .description("Boundary condition type.")
244  .input_selection( get_bc_type_selection() )
245  .input_default("\"none\"")
246  .units( UnitSI::dimensionless() );
247 
248  *this += bc_pressure
249  .disable_where(bc_type, {none, seepage} )
250  .name("bc_pressure")
251  .description("Prescribed pressure value on the boundary. Used for all values of ``bc_type`` except ``none`` and ``seepage``. "
252  "See documentation of ``bc_type`` for exact meaning of ``bc_pressure`` in individual boundary condition types.")
253  .input_default("0.0")
254  .units( UnitSI().m() );
255 
256  *this += bc_flux
257  .disable_where(bc_type, {none, dirichlet} )
258  .name("bc_flux")
259  .description("Incoming water boundary flux. Used for bc_types : ``total_flux``, ``seepage``, ``river``.")
260  .input_default("0.0")
261  .units( UnitSI().m().s(-1) );
262 
263  *this += bc_robin_sigma
264  .disable_where(bc_type, {none, dirichlet, seepage} )
265  .name("bc_robin_sigma")
266  .description("Conductivity coefficient in the ``total_flux`` or the ``river`` boundary condition type.")
267  .input_default("0.0")
268  .units( UnitSI().s(-1) );
269 
270  *this += bc_switch_pressure
271  .disable_where(bc_type, {none, dirichlet, total_flux} )
272  .name("bc_switch_pressure")
273  .description("Critical switch pressure for ``seepage`` and ``river`` boundary conditions.")
274  .input_default("0.0")
275  .units( UnitSI().m() );
276 
277 
278  //these are for unsteady
279  *this += init_pressure.name("init_pressure")
280  .description("Initial condition for pressure in time dependent problems.")
281  .input_default("0.0")
282  .units( UnitSI().m() );
283 
284  *this += storativity.name("storativity")
285  .description("Storativity (in time dependent problems).")
286  .input_default("0.0")
287  .units( UnitSI().m(-1) );
288 
289  *this += extra_storativity.name("extra_storativity")
290  .description("Storativity added from upstream equation.")
291  .units( UnitSI().m(-1) )
292  .input_default("0.0")
293  .flags( input_copy );
294 
295  *this += extra_source.name("extra_water_source_density")
296  .description("Water source density added from upstream equation.")
297  .input_default("0.0")
298  .units( UnitSI().s(-1) )
299  .flags( input_copy );
300 
301  *this += gravity_field.name("gravity")
302  .description("Gravity vector.")
303  .input_default("0.0")
304  .units( UnitSI::dimensionless() );
305 
306  *this += bc_gravity.name("bc_gravity")
307  .description("Boundary gravity vector.")
308  .input_default("0.0")
309  .units( UnitSI::dimensionless() );
310 
311  *this += init_piezo_head.name("init_piezo_head")
312  .units(UnitSI().m())
313  .input_default("0.0")
314  .description("Init piezo head.");
315 
316  *this += bc_piezo_head.name("bc_piezo_head")
317  .units(UnitSI().m())
318  .input_default("0.0")
319  .description("Boundary piezo head.");
320 
321  *this += bc_switch_piezo_head.name("bc_switch_piezo_head")
322  .units(UnitSI().m())
323  .input_default("0.0")
324  .description("Boundary switch piezo head.");
325 
326  //time_term_fields = this->subset({"storativity"});
327  //main_matrix_fields = this->subset({"anisotropy", "conductivity", "cross_section", "sigma", "bc_type", "bc_robin_sigma"});
328  //rhs_fields = this->subset({"water_source_density", "bc_pressure", "bc_flux"});
329 }
330 
332 {
333  mortar_method_=NoMortar;
334 }
335 
336 
337 
338 
339 //=============================================================================
340 // CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL
341 // - do it in parallel:
342 // - initial distribution of elements, edges
343 //
344 /*! @brief CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL
345  *
346  * Parameters {Solver,NSchurs} number of performed Schur
347  * complements (0,1,2) for water flow MH-system
348  *
349  */
350 //=============================================================================
351 DarcyMH::DarcyMH(Mesh &mesh_in, const Input::Record in_rec, TimeGovernor *tm)
352 : DarcyFlowInterface(mesh_in, in_rec),
353  output_object(nullptr),
354  data_changed_(false),
355  schur0(nullptr),
356  steady_diagonal(nullptr),
357  steady_rhs(nullptr),
358  new_diagonal(nullptr),
359  previous_solution(nullptr)
360 {
361 
362  START_TIMER("Darcy constructor");
363  {
364  auto time_rec = in_rec.val<Input::Record>("time");
365  if (tm == nullptr)
366  {
367  time_ = new TimeGovernor(time_rec);
368  }
369  else
370  {
371  TimeGovernor tm_from_rec(time_rec);
372  if (!tm_from_rec.is_default()) // is_default() == false when time record is present in input file
373  {
374  MessageOut() << "Duplicate key 'time', time in flow equation is already initialized from parent class!";
375  ASSERT_PERMANENT(false);
376  }
377  time_ = tm;
378  }
379  }
380 
381  eq_fields_ = make_shared<EqFields>();
382  eq_data_ = make_shared<EqData>();
383  this->eq_fieldset_ = eq_fields_.get();
384 
385  eq_data_->is_linear=true;
386 
387  size = mesh_->n_elements() + mesh_->n_sides() + mesh_->n_edges();
388  n_schur_compls = in_rec.val<int>("n_schurs");
389  eq_data_->mortar_method_= in_rec.val<MortarMethod>("mortar_method");
390  if (eq_data_->mortar_method_ != NoMortar) {
392  }
393 
394 
395 
396  //side_ds->view( std::cout );
397  //el_ds->view( std::cout );
398  //edge_ds->view( std::cout );
399  //rows_ds->view( std::cout );
400 
401 }
402 
404 {
405  // DebugOut() << "t = " << time_->t() << " step_end " << time_->step().end() << "\n";
407  {
408  // In steady case, the solution is computed with the data present at time t,
409  // and the steady state solution is valid until another change in data,
410  // which should correspond to time (t+dt).
411  // "The data change appears immediatly."
412  double next_t = time_->t() + time_->estimate_dt();
413  // DebugOut() << "STEADY next_t = " << next_t << "\n";
414  return next_t * (1 - 2*std::numeric_limits<double>::epsilon());
415  }
416  else
417  {
418  // In unsteady case, the solution is computed with the data present at time t,
419  // and the solution is valid at the time t+dt.
420  // "The data change does not appear immediatly, it is integrated over time interval dt."
421  // DebugOut() << "UNSTEADY\n";
422  return time_->t();
423  }
424 }
425 
427 //connecting data fields with mesh
428 {
429 
430  START_TIMER("data init");
431  eq_data_->mesh = mesh_;
432  eq_fields_->set_mesh(*mesh_);
433 
434  auto gravity_array = input_record_.val<Input::Array>("gravity");
435  std::vector<double> gvec;
436  gravity_array.copy_to(gvec);
437  gvec.push_back(0.0); // zero pressure shift
438  eq_data_->gravity_ = arma::vec(gvec);
439  eq_data_->gravity_vec_ = eq_data_->gravity_.subvec(0,2);
440 
441  FieldValue<3>::VectorFixed gvalue(eq_data_->gravity_vec_);
442  auto field_algo=std::make_shared<FieldConstant<3, FieldValue<3>::VectorFixed>>();
443  field_algo->set_value(gvalue);
444  eq_fields_->gravity_field.set(field_algo, 0.0);
445  eq_fields_->bc_gravity.set(field_algo, 0.0);
446 
447  eq_fields_->bc_pressure.add_factory(
448  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
449  (eq_fields_->bc_gravity, eq_fields_->X(), eq_fields_->bc_piezo_head) );
450  eq_fields_->bc_switch_pressure.add_factory(
451  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
452  (eq_fields_->bc_gravity, eq_fields_->X(), eq_fields_->bc_switch_piezo_head) );
453  eq_fields_->init_pressure.add_factory(
454  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
455  (eq_fields_->gravity_field, eq_fields_->X(), eq_fields_->init_piezo_head) );
456 
457 
458  eq_fields_->set_input_list( this->input_record_.val<Input::Array>("input_fields"), *time_ );
459  // Check that the time step was set for the transient simulation.
460  if (! zero_time_term(true) && time_->is_default() ) {
461  //THROW(ExcAssertMsg());
462  //THROW(ExcMissingTimeGovernor() << input_record_.ei_address());
463  MessageOut() << "Missing the key 'time', obligatory for the transient problems." << endl;
464  ASSERT_PERMANENT(false);
465  }
466 
467  eq_fields_->mark_input_times(*time_);
468 }
469 
470 
471 
473 
474  { // init DOF handler for pressure fields
475 // std::shared_ptr< FiniteElement<0> > fe0_rt = std::make_shared<FE_RT0_disc<0>>();
476  std::shared_ptr< FiniteElement<1> > fe1_rt = std::make_shared<FE_RT0_disc<1>>();
477  std::shared_ptr< FiniteElement<2> > fe2_rt = std::make_shared<FE_RT0_disc<2>>();
478  std::shared_ptr< FiniteElement<3> > fe3_rt = std::make_shared<FE_RT0_disc<3>>();
479  std::shared_ptr< FiniteElement<0> > fe0_disc = std::make_shared<FE_P_disc<0>>(0);
480  std::shared_ptr< FiniteElement<1> > fe1_disc = std::make_shared<FE_P_disc<1>>(0);
481  std::shared_ptr< FiniteElement<2> > fe2_disc = std::make_shared<FE_P_disc<2>>(0);
482  std::shared_ptr< FiniteElement<3> > fe3_disc = std::make_shared<FE_P_disc<3>>(0);
483  std::shared_ptr< FiniteElement<0> > fe0_cr = std::make_shared<FE_CR<0>>();
484  std::shared_ptr< FiniteElement<1> > fe1_cr = std::make_shared<FE_CR<1>>();
485  std::shared_ptr< FiniteElement<2> > fe2_cr = std::make_shared<FE_CR<2>>();
486  std::shared_ptr< FiniteElement<3> > fe3_cr = std::make_shared<FE_CR<3>>();
487 // static FiniteElement<0> fe0_sys = FE_P_disc<0>(0); //TODO fix and use solution with FESystem<0>( {fe0_rt, fe0_disc, fe0_cr} )
488  FESystem<0> fe0_sys( {fe0_disc, fe0_disc, fe0_cr} );
489  FESystem<1> fe1_sys( {fe1_rt, fe1_disc, fe1_cr} );
490  FESystem<2> fe2_sys( {fe2_rt, fe2_disc, fe2_cr} );
491  FESystem<3> fe3_sys( {fe3_rt, fe3_disc, fe3_cr} );
492  MixedPtr<FESystem> fe_sys( std::make_shared<FESystem<0>>(fe0_sys), std::make_shared<FESystem<1>>(fe1_sys),
493  std::make_shared<FESystem<2>>(fe2_sys), std::make_shared<FESystem<3>>(fe3_sys) );
494  std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( mesh_, fe_sys);
495  eq_data_->dh_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
496  eq_data_->dh_->distribute_dofs(ds);
497  }
498 
499  init_eq_data();
500  this->eq_data_->multidim_assembler = AssemblyFlowBase::create< AssemblyMH >(eq_fields_, eq_data_);
502 
503  eq_fields_->add_coords_field();
504 
505  { // construct pressure, velocity and piezo head fields
506  uint rt_component = 0;
507  eq_data_->full_solution = eq_data_->dh_->create_vector();
508  auto ele_flux_ptr = create_field_fe<3, FieldValue<3>::VectorFixed>(eq_data_->dh_, &eq_data_->full_solution, rt_component);
509  eq_fields_->flux.set(ele_flux_ptr, 0.0);
510 
511  eq_fields_->field_ele_velocity.set(Model<3, FieldValue<3>::VectorFixed>::create(fn_mh_velocity(), eq_fields_->flux, eq_fields_->cross_section), 0.0);
512 
513  uint p_ele_component = 1;
514  auto ele_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_, &eq_data_->full_solution, p_ele_component);
515  eq_fields_->field_ele_pressure.set(ele_pressure_ptr, 0.0);
516 
517  uint p_edge_component = 2;
518  auto edge_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_, &eq_data_->full_solution, p_edge_component);
519  eq_fields_->field_edge_pressure.set(edge_pressure_ptr, 0.0);
520 
521  eq_fields_->field_ele_piezo_head.set(
522  Model<3, FieldValue<3>::Scalar>::create(fn_mh_piezohead(), eq_fields_->gravity_field, eq_fields_->X(), eq_fields_->field_ele_pressure),
523  0.0
524  );
525  }
526 
527  { // init DOF handlers represents edge DOFs
528  uint p_edge_component = 2;
529  eq_data_->dh_cr_ = std::make_shared<SubDOFHandlerMultiDim>(eq_data_->dh_,p_edge_component);
530  }
531 
532  { // init DOF handlers represents side DOFs
533  MixedPtr<FE_CR_disc> fe_cr_disc;
534  std::shared_ptr<DiscreteSpace> ds_cr_disc = std::make_shared<EqualOrderDiscreteSpace>( mesh_, fe_cr_disc);
535  eq_data_->dh_cr_disc_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
536  eq_data_->dh_cr_disc_->distribute_dofs(ds_cr_disc);
537  }
538 
539  // Initialize bc_switch_dirichlet to size of global boundary.
540  eq_data_->bc_switch_dirichlet.resize(mesh_->n_elements()+mesh_->bc_mesh()->n_elements(), 1);
541 
542 
545  .val<Input::Record>("nonlinear_solver")
546  .val<Input::AbstractRecord>("linear_solver");
547 
548  // auxiliary set_time call since allocation assembly evaluates fields as well
549  time_->step().use_fparser_ = true;
552 
553 
554 
555  // allocate time term vectors
556  VecDuplicate(schur0->get_solution(), &previous_solution);
557  VecCreateMPI(PETSC_COMM_WORLD, eq_data_->dh_->distr()->lsize(),PETSC_DETERMINE,&(steady_diagonal));
558  VecDuplicate(steady_diagonal,& new_diagonal);
559  VecZeroEntries(new_diagonal);
560  VecDuplicate(steady_diagonal, &steady_rhs);
561 
562 
563  // initialization of balance object
564  balance_ = std::make_shared<Balance>("water", mesh_);
565  balance_->init_from_input(input_record_.val<Input::Record>("balance"), time());
566  eq_data_->water_balance_idx = balance_->add_quantity("water_volume");
567  balance_->allocate(eq_data_->dh_->distr()->lsize(), 1);
568  balance_->units(UnitSI().m(3));
569 
570 
571  eq_data_->balance = balance_;
572  eq_data_->lin_sys = schur0;
573 
574 
576 }
577 
579 {}
580 
582 {
583 
584  /* TODO:
585  * - Allow solution reconstruction (pressure and velocity) from initial condition on user request.
586  * - Steady solution as an intitial condition may be forced by setting inti_time =-1, and set data for the steady solver in that time.
587  * Solver should be able to switch from and to steady case depending on the zero time term.
588  */
589 
590  time_->step().use_fparser_ = true;
592 
593  // zero_time_term means steady case
594  bool zero_time_term_from_right = zero_time_term();
595 
596 
597  if (zero_time_term_from_right) {
598  MessageOut() << "Flow zero time step - steady case\n";
599  // steady case
600  VecZeroEntries(schur0->get_solution());
601  //read_initial_condition(); // Possible solution guess for steady case.
602  use_steady_assembly_ = true;
603  solve_nonlinear(); // with right limit data
604  // eq_data_->full_solution.local_to_ghost_begin();
605  // eq_data_->full_solution.local_to_ghost_end();
606  } else {
607  MessageOut() << "Flow zero time step - unsteady case\n";
608  VecZeroEntries(schur0->get_solution());
609  VecZeroEntries(previous_solution);
610 
612  reconstruct_solution_from_schur(eq_data_->multidim_assembler);
613 
614  assembly_linear_system(); // in particular due to balance
615 // print_matlab_matrix("matrix_zero");
616  }
617  //solution_output(T,right_limit); // data for time T in any case
618 
619  eq_data_->full_solution.local_to_ghost_begin();
620  eq_data_->full_solution.local_to_ghost_end();
621  output_data();
622 }
623 
624 //=============================================================================
625 // COMPOSE and SOLVE WATER MH System possibly through Schur complements
626 //=============================================================================
628 {
629  START_TIMER("Solving MH system");
630 
631  time_->next_time();
632 
633  time_->view("DARCY"); //time governor information output
634 
635  solve_time_step();
636 
637  eq_data_->full_solution.local_to_ghost_begin();
638  eq_data_->full_solution.local_to_ghost_end();
639 }
640 
641 
642 void DarcyMH::solve_time_step(bool output)
643 {
644  time_->step().use_fparser_ = true;
646 
647  bool zero_time_term_from_left=zero_time_term();
648 
649  bool jump_time = eq_fields_->storativity.is_jump_time();
650  if (! zero_time_term_from_left) {
651  MessageOut() << "Flow time step - unsteady case\n";
652  // time term not treated as zero
653  // Unsteady solution up to the T.
654 
655  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
656  use_steady_assembly_ = false;
657 
659  solve_nonlinear(); // with left limit data
660  if (jump_time) {
661  WarningOut() << "Output of solution discontinuous in time not supported yet.\n";
662  //solution_output(T, left_limit); // output use time T- delta*dt
663  //output_data();
664  }
665  }
666 
667  if (time_->is_end()) {
668  // output for unsteady case, end_time should not be the jump time
669  // but rether check that
670  if (! zero_time_term_from_left && ! jump_time && output) output_data();
671  return;
672  }
673 
674  time_->step().use_fparser_ = true;
676  bool zero_time_term_from_right=zero_time_term();
677  if (zero_time_term_from_right) {
678  MessageOut() << "Flow time step - steady case\n";
679  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
680  use_steady_assembly_ = true;
681  solve_nonlinear(); // with right limit data
682 
683  } else if (! zero_time_term_from_left && jump_time) {
684  WarningOut() << "Discontinuous time term not supported yet.\n";
685  //solution_transfer(); // internally call set_time(T, left) and set_time(T,right) again
686  //solve_nonlinear(); // with right limit data
687  }
688 
689  //solution_output(T,right_limit); // data for time T in any case
690  if (output) output_data();
691 }
692 
693 
694 bool DarcyMH::zero_time_term(bool time_global) {
695  if (time_global) {
696  return (eq_fields_->storativity.input_list_size() == 0);
697  } else {
698  return eq_fields_->storativity.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
699  }
700 }
701 
702 
704 {
705 
707  double residual_norm = schur0->compute_residual();
709  MessageOut().fmt("[nonlinear solver] norm of initial residual: {}\n", residual_norm);
710 
711  // Reduce is_linear flag.
712  int is_linear_common;
713  MPI_Allreduce(&(eq_data_->is_linear), &is_linear_common,1, MPI_INT ,MPI_MIN,PETSC_COMM_WORLD);
714 
715  Input::Record nl_solver_rec = input_record_.val<Input::Record>("nonlinear_solver");
716  this->tolerance_ = nl_solver_rec.val<double>("tolerance");
717  this->max_n_it_ = nl_solver_rec.val<unsigned int>("max_it");
718  this->min_n_it_ = nl_solver_rec.val<unsigned int>("min_it");
719  if (this->min_n_it_ > this->max_n_it_) this->min_n_it_ = this->max_n_it_;
720 
721  if (! is_linear_common) {
722  // set tolerances of the linear solver unless they are set by user.
723  schur0->set_tolerances(0.1*this->tolerance_, 0.01*this->tolerance_, 100);
724  }
725  vector<double> convergence_history;
726 
727  Vec save_solution;
728  VecDuplicate(schur0->get_solution(), &save_solution);
729  while (nonlinear_iteration_ < this->min_n_it_ ||
730  (residual_norm > this->tolerance_ && nonlinear_iteration_ < this->max_n_it_ )) {
731  ASSERT_EQ( convergence_history.size(), nonlinear_iteration_ );
732  convergence_history.push_back(residual_norm);
733 
734  // stagnation test
735  if (convergence_history.size() >= 5 &&
736  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 2] > 0.9 &&
737  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 5] > 0.8) {
738  // stagnation
739  if (input_record_.val<Input::Record>("nonlinear_solver").val<bool>("converge_on_stagnation")) {
740  WarningOut().fmt("Accept solution on stagnation. Its: {} Residual: {}\n", nonlinear_iteration_, residual_norm);
741  break;
742  } else {
743  THROW(ExcSolverDiverge() << EI_Reason("Stagnation."));
744  }
745  }
746 
747  if (! is_linear_common)
748  VecCopy( schur0->get_solution(), save_solution);
751 
752  // hack to make BDDC work with empty compute_residual
753  if (is_linear_common){
754  // we want to print this info in linear (and steady) case
755  residual_norm = schur0->compute_residual();
756  MessageOut().fmt("[nonlinear solver] lin. it: {}, reason: {}, residual: {}\n",
757  si.n_iterations, si.converged_reason, residual_norm);
758  break;
759  }
760  data_changed_=true; // force reassembly for non-linear case
761 
762  double alpha = 1; // how much of new solution
763  VecAXPBY(schur0->get_solution(), (1-alpha), alpha, save_solution);
764 
765  //LogOut().fmt("Linear solver ended with reason: {} \n", si.converged_reason );
766  //ASSERT_PERMANENT_GE( si.converged_reason, 0).error("Linear solver failed to converge. \n");
768 
769  residual_norm = schur0->compute_residual();
770  MessageOut().fmt("[nonlinear solver] it: {} lin. it: {}, reason: {}, residual: {}\n",
771  nonlinear_iteration_, si.n_iterations, si.converged_reason, residual_norm);
772  }
773  chkerr(VecDestroy(&save_solution));
774  this -> postprocess();
775 
776  // adapt timestep
777  if (! this->zero_time_term()) {
778  double mult = 1.0;
779  if (nonlinear_iteration_ < 3) mult = 1.6;
780  if (nonlinear_iteration_ > 7) mult = 0.7;
781  time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
782  // int result = time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
783  //DebugOut().fmt("time adaptivity, res: {} it: {} m: {} dt: {} edt: {}\n", result, nonlinear_iteration_, mult, time_->dt(), time_->estimate_dt());
784  }
785 }
786 
788 {
790 }
791 
793 {
794  START_TIMER("postprocess");
795 
796  //fix velocity when mortar method is used
797  if(eq_data_->mortar_method_ != MortarMethod::NoMortar){
798  auto multidim_assembler = AssemblyFlowBase::create< AssemblyMH >(eq_fields_, eq_data_);
799  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
800  unsigned int dim = dh_cell.dim();
801  multidim_assembler[dim-1]->fix_velocity(dh_cell);
802  }
803  }
804  //ElementAccessor<3> ele;
805 
806  // modify side fluxes in parallel
807  // for every local edge take time term on digonal and add it to the corresponding flux
808  /*
809  for (unsigned int i_loc = 0; i_loc < el_ds->lsize(); i_loc++) {
810  ele = mesh_->element_accessor(el_4_loc[i_loc]);
811  for (unsigned int i=0; i<ele->n_sides(); i++) {
812  side_rows[i] = side_row_4_id[ mh_dh.side_dof( ele_ac.side(i) ) ];
813  values[i] = -1.0 * ele_ac.measure() *
814  data.cross_section.value(ele_ac.centre(), ele_ac.element_accessor()) *
815  data.water_source_density.value(ele_ac.centre(), ele_ac.element_accessor()) /
816  ele_ac.n_sides();
817  }
818  VecSetValues(schur0->get_solution(), ele_ac.n_sides(), side_rows, values, ADD_VALUES);
819  }
820  VecAssemblyBegin(schur0->get_solution());
821  VecAssemblyEnd(schur0->get_solution());
822  */
823 }
824 
825 
827  START_TIMER("Darcy output data");
828 
829  //print_matlab_matrix("matrix_" + std::to_string(time_->step().index()));
830 
831  //time_->view("DARCY"); //time governor information output
832  this->output_object->output();
833 
834 
835  START_TIMER("Darcy balance output");
836  balance_->calculate_cumulative(eq_data_->water_balance_idx, schur0->get_solution());
837  balance_->calculate_instant(eq_data_->water_balance_idx, schur0->get_solution());
838  balance_->output();
839 }
840 
841 
843 {
844  return schur0->get_solution_precision();
845 }
846 
847 
848 // ===========================================================================================
849 //
850 // MATRIX ASSEMBLY - we use abstract assembly routine, where LS Mat/Vec SetValues
851 // are in fact pointers to allocating or filling functions - this is governed by Linsystem roitunes
852 //
853 // =======================================================================================
855 {
856  START_TIMER("DarcyFlowMHy::assembly_mh_matrix");
857 
858  // set auxiliary flag for switchting Dirichlet like BC
859  eq_data_->force_no_neumann_bc = use_steady_assembly_ && (nonlinear_iteration_ == 0);
860  eq_data_->n_schur_compls = n_schur_compls;
861 
862 
863  balance_->start_flux_assembly(eq_data_->water_balance_idx);
864 
865  // TODO: try to move this into balance, or have it in the generic assembler class, that should perform the cell loop
866  // including various pre- and post-actions
867  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
868  unsigned int dim = dh_cell.dim();
869  assembler[dim-1]->assemble(dh_cell);
870  }
871 
872 
873  balance_->finish_flux_assembly(eq_data_->water_balance_idx);
874 
875 }
876 
877 
879 {
880  START_TIMER("DarcyFlowMH::allocate_mh_matrix");
881 
882  // set auxiliary flag for switchting Dirichlet like BC
883  eq_data_->n_schur_compls = n_schur_compls;
884  LinSys *ls = schur0;
885 
886  // to make space for second schur complement, max. 10 neighbour edges of one el.
887  double zeros[100000];
888  for(int i=0; i<100000; i++) zeros[i] = 0.0;
889 
890  std::vector<LongIdx> tmp_rows;
891  tmp_rows.reserve(200);
892 
893  std::vector<LongIdx> dofs, dofs_ngh;
894  dofs.reserve(eq_data_->dh_->max_elem_dofs());
895  dofs_ngh.reserve(eq_data_->dh_->max_elem_dofs());
896 
897  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
898  ElementAccessor<3> ele = dh_cell.elm();
899 
900  const uint ndofs = dh_cell.n_dofs();
901  dofs.resize(ndofs);
902  dh_cell.get_dof_indices(dofs);
903 
904  // whole local MH matrix
905  ls->mat_set_values(ndofs, dofs.data(), ndofs, dofs.data(), zeros);
906 
907  tmp_rows.clear();
908 
909  // compatible neighborings rows
910  unsigned int n_neighs = ele->n_neighs_vb();
911  for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
912  // every compatible connection adds a 2x2 matrix involving
913  // current element pressure and a connected edge pressure
914 
915  // read neighbor dofs (dh dofhandler)
916  // neighbor cell owning neighb_side
917  DHCellAccessor dh_neighb_cell = neighb_side.cell();
918 
919  const uint ndofs_ngh = dh_neighb_cell.n_dofs();
920  dofs_ngh.resize(ndofs_ngh);
921  dh_neighb_cell.get_dof_indices(dofs_ngh);
922 
923  // local index of pedge dof on neighboring cell
924  // (dim+1) is number of edges of higher dim element
925  // TODO: replace with DHCell getter when available for FESystem component
926  const unsigned int t = dh_neighb_cell.n_dofs() - (dh_neighb_cell.dim()+1) + neighb_side.side().side_idx();
927  tmp_rows.push_back(dofs_ngh[t]);
928  }
929 
930  const uint nsides = ele->n_sides();
931  LongIdx * edge_rows = dofs.data() + nsides; // pointer to start of ele
932  // allocate always also for schur 2
933  ls->mat_set_values(nsides+1, edge_rows, n_neighs, tmp_rows.data(), zeros); // (edges, ele) x (neigh edges)
934  ls->mat_set_values(n_neighs, tmp_rows.data(), nsides+1, edge_rows, zeros); // (neigh edges) x (edges, ele)
935  ls->mat_set_values(n_neighs, tmp_rows.data(), n_neighs, tmp_rows.data(), zeros); // (neigh edges) x (neigh edges)
936 
937  tmp_rows.clear();
938 
939  if (eq_data_->mortar_method_ != NoMortar) {
940  auto &isec_list = mesh_->mixed_intersections().element_intersections_[ele.idx()];
941  for(auto &isec : isec_list ) {
942  IntersectionLocalBase *local = isec.second;
943  DHCellAccessor dh_cell_slave = eq_data_->dh_->cell_accessor_from_element(local->bulk_ele_idx());
944 
945  const uint ndofs_slave = dh_cell_slave.n_dofs();
946  dofs_ngh.resize(ndofs_slave);
947  dh_cell_slave.get_dof_indices(dofs_ngh);
948 
949  //DebugOut().fmt("Alloc: {} {}", ele.idx(), local->bulk_ele_idx());
950  for(unsigned int i_side=0; i_side < dh_cell_slave.elm()->n_sides(); i_side++) {
951  // TODO: replace with DHCell getter when available for FESystem component
952  tmp_rows.push_back( dofs_ngh[(ndofs_slave+1)/2+i_side] );
953  //DebugOut() << "aedge" << print_var(tmp_rows[tmp_rows.size()-1]);
954  }
955  }
956  }
957  /*
958  for(unsigned int i_side=0; i_side < ele->n_sides(); i_side++) {
959  DebugOut() << "aedge:" << print_var(edge_rows[i_side]);
960  }*/
961 
962  edge_rows = dofs.data() + nsides +1; // pointer to start of edges
963  ls->mat_set_values(nsides, edge_rows, tmp_rows.size(), tmp_rows.data(), zeros); // master edges x neigh edges
964  ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), nsides, edge_rows, zeros); // neigh edges x master edges
965  ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), tmp_rows.size(), tmp_rows.data(), zeros); // neigh edges x neigh edges
966 
967  }
968 /*
969  // alloc edge diagonal entries
970  if(rank == 0)
971  for( vector<Edge>::iterator edg = mesh_->edges.begin(); edg != mesh_->edges.end(); ++edg) {
972  int edg_idx = mh_dh.row_4_edge[edg->side(0)->edge_idx()];
973 
974 // for( vector<Edge>::iterator edg2 = mesh_->edges.begin(); edg2 != mesh_->edges.end(); ++edg2){
975 // int edg_idx2 = mh_dh.row_4_edge[edg2->side(0)->edge_idx()];
976 // if(edg_idx == edg_idx2){
977 // DBGCOUT(<< "P[ " << rank << " ] " << "edg alloc: " << edg_idx << " " << edg_idx2 << "\n");
978  ls->mat_set_value(edg_idx, edg_idx, 0.0);
979 // }
980 // }
981  }
982  */
983  /*
984  if (mortar_method_ == MortarP0) {
985  P0_CouplingAssembler(*this).assembly(*ls);
986  } else if (mortar_method_ == MortarP1) {
987  P1_CouplingAssembler(*this).assembly(*ls);
988  }*/
989 }
990 
992 {
993  START_TIMER("assembly source term");
994  balance_->start_source_assembly(eq_data_->water_balance_idx);
995 
996  std::vector<LongIdx> global_dofs(eq_data_->dh_->max_elem_dofs());
997 
998  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
999  ElementAccessor<3> ele = dh_cell.elm();
1000 
1001  const uint ndofs = dh_cell.n_dofs();
1002  global_dofs.resize(ndofs);
1003  dh_cell.get_dof_indices(global_dofs);
1004 
1005  double cs = eq_fields_->cross_section.value(ele.centre(), ele);
1006 
1007  // set sources
1008  double source = ele.measure() * cs *
1009  (eq_fields_->water_source_density.value(ele.centre(), ele)
1010  +eq_fields_->extra_source.value(ele.centre(), ele));
1011  // TODO: replace with DHCell getter when available for FESystem component
1012  schur0->rhs_set_value(global_dofs[ndofs/2], -1.0 * source );
1013 
1014  balance_->add_source_values(eq_data_->water_balance_idx, ele.region().bulk_idx(),
1015  {dh_cell.get_loc_dof_indices()[ndofs/2]}, {0}, {source});
1016  }
1017 
1018  balance_->finish_source_assembly(eq_data_->water_balance_idx);
1019 }
1020 
1021 
1022 
1023 
1024 /*******************************************************************************
1025  * COMPOSE WATER MH MATRIX WITHOUT SCHUR COMPLEMENT
1026  ******************************************************************************/
1027 
1029 
1030  START_TIMER("preallocation");
1031 
1032  if (schur0 == NULL) { // create Linear System for MH matrix
1033 
1034  if (in_rec.type() == LinSys_BDDC::get_input_type()) {
1035 #ifdef FLOW123D_HAVE_BDDCML
1036  WarningOut() << "For BDDC no Schur complements are used.";
1037  n_schur_compls = 0;
1038  LinSys_BDDC *ls = new LinSys_BDDC(&(*eq_data_->dh_->distr()),
1039  true); // swap signs of matrix and rhs to make the matrix SPD
1040  ls->set_from_input(in_rec);
1041  ls->set_solution( eq_data_->full_solution.petsc_vec() );
1042  // possible initialization particular to BDDC
1043  START_TIMER("BDDC set mesh data");
1045  schur0=ls;
1046  END_TIMER("BDDC set mesh data");
1047 #else
1048  //Exception
1049  THROW( ExcBddcmlNotSupported() );
1050 #endif // FLOW123D_HAVE_BDDCML
1051  }
1052  else if (in_rec.type() == LinSys_PETSC::get_input_type()) {
1053  // use PETSC for serial case even when user wants BDDC
1054  if (n_schur_compls > 2) {
1055  WarningOut() << "Invalid number of Schur Complements. Using 2.";
1056  n_schur_compls = 2;
1057  }
1058 
1059  LinSys_PETSC *schur1, *schur2;
1060 
1061  if (n_schur_compls == 0) {
1062  LinSys_PETSC *ls = new LinSys_PETSC( &(*eq_data_->dh_->distr()) );
1063 
1064  // temporary solution; we have to set precision also for sequantial case of BDDC
1065  // final solution should be probably call of direct solver for oneproc case
1066  if (in_rec.type() != LinSys_BDDC::get_input_type()) ls->set_from_input(in_rec);
1067  else {
1068  ls->LinSys::set_from_input(in_rec); // get only common options
1069  }
1070 
1071  ls->set_solution( eq_data_->full_solution.petsc_vec() );
1072  schur0=ls;
1073  } else {
1074  IS is;
1075  auto side_dofs_vec = get_component_indices_vec(0);
1076 
1077  ISCreateGeneral(PETSC_COMM_SELF, side_dofs_vec.size(), &(side_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1078  //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1079  //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1080 
1081  SchurComplement *ls = new SchurComplement(&(*eq_data_->dh_->distr()), is);
1082 
1083  // make schur1
1085  if (n_schur_compls==1) {
1086  schur1 = new LinSys_PETSC(ds);
1087  schur1->set_positive_definite();
1088  } else {
1089  IS is;
1090  auto elem_dofs_vec = get_component_indices_vec(1);
1091 
1092  const PetscInt *b_indices;
1093  ISGetIndices(ls->IsB, &b_indices);
1094  uint b_size = ls->loc_size_B;
1095  for(uint i_b=0, i_bb=0; i_b < b_size && i_bb < elem_dofs_vec.size(); i_b++) {
1096  if (b_indices[i_b] == elem_dofs_vec[i_bb])
1097  elem_dofs_vec[i_bb++] = i_b + ds->begin();
1098  }
1099  ISRestoreIndices(ls->IsB, &b_indices);
1100 
1101 
1102  ISCreateGeneral(PETSC_COMM_SELF, elem_dofs_vec.size(), &(elem_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1103  //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1104  //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1105  SchurComplement *ls1 = new SchurComplement(ds, is); // is is deallocated by SchurComplement
1106  ls1->set_negative_definite();
1107 
1108  // make schur2
1109  schur2 = new LinSys_PETSC( ls1->make_complement_distribution() );
1110  schur2->set_positive_definite();
1111  ls1->set_complement( schur2 );
1112  schur1 = ls1;
1113  }
1114  ls->set_complement( schur1 );
1115  ls->set_from_input(in_rec);
1116  ls->set_solution( eq_data_->full_solution.petsc_vec() );
1117  schur0=ls;
1118  }
1119 
1120  START_TIMER("PETSC PREALLOCATION");
1121  schur0->set_symmetric();
1123 
1125 
1126  VecZeroEntries(schur0->get_solution());
1127  END_TIMER("PETSC PREALLOCATION");
1128  }
1129  else {
1130  THROW( ExcUnknownSolver() );
1131  }
1132 
1133  END_TIMER("preallocation");
1134  }
1135 
1136 }
1137 
1138 
1139 
1140 
1142  START_TIMER("DarcyFlowMH_Steady::assembly_linear_system");
1143 
1144  eq_data_->is_linear=true;
1145  bool is_steady = zero_time_term();
1146  //DebugOut() << "Assembly linear system\n";
1147  if (data_changed_) {
1148  data_changed_ = false;
1149  //DebugOut() << "Data changed\n";
1150  // currently we have no optimization for cases when just time term data or RHS data are changed
1151  START_TIMER("full assembly");
1152  if (typeid(*schur0) != typeid(LinSys_BDDC)) {
1153  schur0->start_add_assembly(); // finish allocation and create matrix
1154  }
1155 
1158 
1160 
1161 
1162  assembly_mh_matrix( eq_data_->multidim_assembler ); // fill matrix
1163 
1165 // print_matlab_matrix("matrix");
1167  //MatView( *const_cast<Mat*>(schur0->get_matrix()), PETSC_VIEWER_STDOUT_WORLD );
1168  //VecView( *const_cast<Vec*>(schur0->get_rhs()), PETSC_VIEWER_STDOUT_WORLD);
1169 
1170  if (! is_steady) {
1171  START_TIMER("fix time term");
1172  //DebugOut() << "setup time term\n";
1173  // assembly time term and rhs
1174  setup_time_term();
1175  modify_system();
1176  }
1177  else
1178  {
1179  balance_->start_mass_assembly(eq_data_->water_balance_idx);
1180  balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1181  }
1182  END_TIMER("full assembly");
1183  } else {
1184  START_TIMER("modify system");
1185  if (! is_steady) {
1186  modify_system();
1187  } else {
1188  //Should be replaced with exception if error will be switched on.
1189  //ASSERT_PERMANENT(false).error("Planned computation time for steady solver, but data are not changed.\n");
1190  }
1191  END_TIMER("modiffy system");
1192  }
1193 
1194 }
1195 
1196 
1197 void DarcyMH::print_matlab_matrix(std::string matlab_file)
1198 {
1199  std::string output_file;
1200 
1201  if ( typeid(*schur0) == typeid(LinSys_BDDC) ){
1202 // WarningOut() << "Can output matrix only on a single processor.";
1203 // output_file = FilePath(matlab_file + "_bddc.m", FilePath::output_file);
1204 // ofstream os( output_file );
1205 // auto bddc = static_cast<LinSys_BDDC*>(schur0);
1206 // bddc->print_matrix(os);
1207  }
1208  else {//if ( typeid(*schur0) == typeid(LinSys_PETSC) ){
1209  output_file = FilePath(matlab_file + ".m", FilePath::output_file);
1210  PetscViewer viewer;
1211  PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file.c_str(), &viewer);
1212  PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1213  MatView( *const_cast<Mat*>(schur0->get_matrix()), viewer);
1214  VecView( *const_cast<Vec*>(schur0->get_rhs()), viewer);
1215  VecView( *const_cast<Vec*>(schur0->get_rhs()), viewer);
1216  VecView( *const_cast<Vec*>(&(schur0->get_solution())), viewer);
1217  }
1218 // else{
1219 // WarningOut() << "No matrix output available for the current solver.";
1220 // return;
1221 // }
1222 
1223  // compute h_min for different dimensions
1224  double d_max = std::numeric_limits<double>::max();
1225  double h1 = d_max, h2 = d_max, h3 = d_max;
1226  double he2 = d_max, he3 = d_max;
1227  for (auto ele : mesh_->elements_range()) {
1228  switch(ele->dim()){
1229  case 1: h1 = std::min(h1,ele.measure()); break;
1230  case 2: h2 = std::min(h2,ele.measure()); break;
1231  case 3: h3 = std::min(h3,ele.measure()); break;
1232  }
1233 
1234  for (unsigned int j=0; j<ele->n_sides(); j++) {
1235  switch(ele->dim()){
1236  case 2: he2 = std::min(he2, ele.side(j)->measure()); break;
1237  case 3: he3 = std::min(he3, ele.side(j)->measure()); break;
1238  }
1239  }
1240  }
1241  if(h1 == d_max) h1 = 0;
1242  if(h2 == d_max) h2 = 0;
1243  if(h3 == d_max) h3 = 0;
1244  if(he2 == d_max) he2 = 0;
1245  if(he3 == d_max) he3 = 0;
1246 
1247  FILE * file;
1248  file = fopen(output_file.c_str(),"a");
1249  fprintf(file, "nA = %d;\n", eq_data_->dh_cr_disc_->distr()->size());
1250  fprintf(file, "nB = %d;\n", eq_data_->dh_->mesh()->get_el_ds()->size());
1251  fprintf(file, "nBF = %d;\n", eq_data_->dh_cr_->distr()->size());
1252  fprintf(file, "h1 = %e;\nh2 = %e;\nh3 = %e;\n", h1, h2, h3);
1253  fprintf(file, "he2 = %e;\nhe3 = %e;\n", he2, he3);
1254  fclose(file);
1255 }
1256 
1257 
1258 //template <int dim>
1259 //std::vector<arma::vec3> dof_points(DHCellAccessor cell, const Mapping<dim, 3> &mapping) {
1260 //
1261 //
1262 // vector<arma::vec::fixed<dim+1>> bary_dof_points = cell->fe()->dof_points();
1263 //
1264 // std::vector<arma::vec3> points(20);
1265 // points.resize(0);
1266 //
1267 //}
1268 //
1269 
1271  START_TIMER("DarcyFlowMH_Steady::set_mesh_data_for_bddc");
1272  // prepare mesh for BDDCML
1273  // initialize arrays
1274  // auxiliary map for creating coordinates of local dofs and global-to-local numbering
1275  std::map<int, arma::vec3> localDofMap;
1276  // connectivity for the subdomain, i.e. global dof numbers on element, stored element-by-element
1277  // Indices of Nodes on Elements
1278  std::vector<int> inet;
1279  // number of degrees of freedom on elements - determines elementwise chunks of INET array
1280  // Number of Nodes on Elements
1281  std::vector<int> nnet;
1282  // Indices of Subdomain Elements in Global Numbering - for local elements, their global indices
1283  std::vector<int> isegn;
1284 
1285  // This array is currently not used in BDDCML, it was used as an interface scaling alternative to scaling
1286  // by diagonal. It corresponds to the rho-scaling.
1287  std::vector<double> element_permeability;
1288 
1289  // maximal and minimal dimension of elements
1290  uint elDimMax = 1;
1291  uint elDimMin = 3;
1292  std::vector<LongIdx> cell_dofs_global(10);
1293 
1294 
1295 
1296  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1297  // for each element, create local numbering of dofs as fluxes (sides), pressure (element centre), Lagrange multipliers (edges), compatible connections
1298 
1299  dh_cell.get_dof_indices(cell_dofs_global);
1300 
1301  inet.insert(inet.end(), cell_dofs_global.begin(), cell_dofs_global.end());
1302  uint n_inet = cell_dofs_global.size();
1303 
1304 
1305  uint dim = dh_cell.elm().dim();
1306  elDimMax = std::max( elDimMax, dim );
1307  elDimMin = std::min( elDimMin, dim );
1308 
1309  // TODO: this is consistent with previous implementation, but may be wrong as it use global element numbering
1310  // used in sequential mesh, do global numbering of distributed elements.
1311  isegn.push_back( dh_cell.elm_idx());
1312 
1313  // TODO: use FiniteElement::dof_points
1314  for (unsigned int si=0; si<dh_cell.elm()->n_sides(); si++) {
1315  arma::vec3 coord = dh_cell.elm().side(si)->centre();
1316  // TODO: replace with DHCell getter when available for FESystem component
1317  // flux dof points
1318  localDofMap.insert( std::make_pair( cell_dofs_global[si], coord ) );
1319  // pressure trace dof points
1320  localDofMap.insert( std::make_pair( cell_dofs_global[si+dim+2], coord ) );
1321  }
1322  // pressure dof points
1323  arma::vec3 elm_centre = dh_cell.elm().centre();
1324  localDofMap.insert( std::make_pair( cell_dofs_global[dim+1], elm_centre ) );
1325 
1326  // insert dofs related to compatible connections
1327  //const Element *ele = dh_cell.elm().element();
1328  for(DHCellSide side : dh_cell.neighb_sides()) {
1329  uint neigh_dim = side.cell().elm().dim();
1330  side.cell().get_dof_indices(cell_dofs_global);
1331  int edge_row = cell_dofs_global[neigh_dim+2+side.side_idx()];
1332  localDofMap.insert( std::make_pair( edge_row, side.centre() ) );
1333  inet.push_back( edge_row );
1334  n_inet++;
1335  }
1336  nnet.push_back(n_inet);
1337 
1338 
1339  // version for rho scaling
1340  // trace computation
1341  double conduct = eq_fields_->conductivity.value( elm_centre , dh_cell.elm() );
1342  auto aniso = eq_fields_->anisotropy.value( elm_centre , dh_cell.elm() );
1343 
1344  // compute mean on the diagonal
1345  double coef = 0.;
1346  for ( int i = 0; i < 3; i++) {
1347  coef = coef + aniso.at(i,i);
1348  }
1349  // Maybe divide by cs
1350  coef = conduct*coef / 3;
1351 
1352  ASSERT_GT( coef, 0. ).error("Zero coefficient of hydrodynamic resistance. \n");
1353  element_permeability.push_back( 1. / coef );
1354  }
1355 // uint i_inet = 0;
1356 // for(int n_dofs : nnet) {
1357 // DebugOut() << "nnet: " << n_dofs;
1358 // for(int j=0; j < n_dofs; j++, i_inet++) {
1359 // DebugOut() << "inet: " << inet[i_inet];
1360 // }
1361 // }
1362 
1363  auto distr = eq_data_->dh_->distr();
1364 // for(auto pair : localDofMap) {
1365 // DebugOut().every_proc() << "r: " << distr->myp() << " gi: " << pair.first << "xyz: " << pair.second[0];
1366 //
1367 // }
1368 
1369 
1370  //convert set of dofs to vectors
1371  // number of nodes (= dofs) on the subdomain
1372  int numNodeSub = localDofMap.size();
1373  //ASSERT_PERMANENT_EQ( (unsigned int)numNodeSub, eq:data_->dh_->lsize() );
1374  // Indices of Subdomain Nodes in Global Numbering - for local nodes, their global indices
1375  std::vector<int> isngn( numNodeSub );
1376  // pseudo-coordinates of local nodes (i.e. dofs)
1377  // they need not be exact, they are used just for some geometrical considerations in BDDCML,
1378  // such as selection of corners maximizing area of a triangle, bounding boxes fro subdomains to
1379  // find candidate neighbours etc.
1380  std::vector<double> xyz( numNodeSub * 3 ) ;
1381  int ind = 0;
1382  std::map<int,arma::vec3>::iterator itB = localDofMap.begin();
1383  for ( ; itB != localDofMap.end(); ++itB ) {
1384  isngn[ind] = itB -> first;
1385 
1386  arma::vec3 coord = itB -> second;
1387  for ( int j = 0; j < 3; j++ ) {
1388  xyz[ j*numNodeSub + ind ] = coord[j];
1389  }
1390 
1391  ind++;
1392  }
1393  localDofMap.clear();
1394 
1395  // Number of Nodal Degrees of Freedom
1396  // nndf is trivially one - dofs coincide with nodes
1397  std::vector<int> nndf( numNodeSub, 1 );
1398 
1399  // prepare auxiliary map for renumbering nodes
1400  typedef std::map<int,int> Global2LocalMap_; //! type for storage of global to local map
1401  Global2LocalMap_ global2LocalNodeMap;
1402  for ( unsigned ind = 0; ind < isngn.size(); ++ind ) {
1403  global2LocalNodeMap.insert( std::make_pair( static_cast<unsigned>( isngn[ind] ), ind ) );
1404  }
1405 
1406  // renumber nodes in the inet array to locals
1407  int indInet = 0;
1408  for ( unsigned int iEle = 0; iEle < isegn.size(); iEle++ ) {
1409  int nne = nnet[ iEle ];
1410  for ( int ien = 0; ien < nne; ien++ ) {
1411 
1412  int indGlob = inet[indInet];
1413  // map it to local node
1414  Global2LocalMap_::iterator pos = global2LocalNodeMap.find( indGlob );
1415  ASSERT( pos != global2LocalNodeMap.end() )(indGlob).error("Cannot remap node global index to local indices.\n");
1416  int indLoc = static_cast<int> ( pos -> second );
1417 
1418  // store the node
1419  inet[ indInet++ ] = indLoc;
1420  }
1421  }
1422 
1423  int numNodes = size;
1424  int numDofsInt = size;
1425  int spaceDim = 3; // TODO: what is the proper value here?
1426  int meshDim = elDimMax;
1427 
1428  /**
1429  * We need:
1430  * - local to global element map (possibly mesh->el_4_loc
1431  * - inet, nnet - local dof numbers per element, local numbering of only those dofs that are on owned elements
1432  * 1. collect DH local dof indices on elements, manage map from DH local indices to BDDC local dof indices
1433  * 2. map collected DH indices to BDDC indices using the map
1434  * - local BDDC dofs to global dofs, use DH to BDDC map with DH local to global map
1435  * - XYZ - permuted, collect in main loop into array of size of all DH local dofs, compress and rearrange latter
1436  * - element_permeability - in main loop
1437  */
1438  bddc_ls -> load_mesh( LinSys_BDDC::BDDCMatrixType::SPD_VIA_SYMMETRICGENERAL, spaceDim, numNodes, numDofsInt, inet, nnet, nndf, isegn, isngn, isngn, xyz, element_permeability, meshDim );
1439 }
1440 
1441 
1442 
1443 
1444 //=============================================================================
1445 // DESTROY WATER MH SYSTEM STRUCTURE
1446 //=============================================================================
1448  if (previous_solution != nullptr) chkerr(VecDestroy(&previous_solution));
1449  if (steady_diagonal != nullptr) chkerr(VecDestroy(&steady_diagonal));
1450  if (new_diagonal != nullptr) chkerr(VecDestroy(&new_diagonal));
1451  if (steady_rhs != nullptr) chkerr(VecDestroy(&steady_rhs));
1452 
1453 
1454  if (schur0 != NULL) {
1455  delete schur0;
1456  }
1457 
1458  if (output_object) delete output_object;
1459 
1460  if(time_ != nullptr)
1461  delete time_;
1462 
1463 }
1464 
1465 
1466 /*
1467 void mat_count_off_proc_values(Mat m, Vec v) {
1468  int n, first, last;
1469  const PetscInt *cols;
1470  Distribution distr(v);
1471 
1472  int n_off = 0;
1473  int n_on = 0;
1474  int n_off_rows = 0;
1475  MatGetOwnershipRange(m, &first, &last);
1476  for (int row = first; row < last; row++) {
1477  MatGetRow(m, row, &n, &cols, PETSC_NULL);
1478  bool exists_off = false;
1479  for (int i = 0; i < n; i++)
1480  if (distr.get_proc(cols[i]) != distr.myp())
1481  n_off++, exists_off = true;
1482  else
1483  n_on++;
1484  if (exists_off)
1485  n_off_rows++;
1486  MatRestoreRow(m, row, &n, &cols, PETSC_NULL);
1487  }
1488 }
1489 */
1490 
1491 
1492 
1493 
1494 
1495 
1496 
1497 
1498 
1499 
1500 
1501 
1503 {
1504  DebugOut().fmt("Setup with dt: {}\n", time_->dt());
1505  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1506  ElementAccessor<3> ele = dh_cell.elm();
1507 
1508  double init_value = eq_fields_->init_pressure.value(ele.centre(),ele);
1509 
1510  LocDofVec loc_dofs = dh_cell.get_loc_dof_indices();
1511  for (unsigned int i=0; i<ele->n_sides(); i++) {
1512  uint n_sides_of_edge = ele.side(i)->edge().n_sides();
1513  uint edge_dof = ele->n_sides()+1+i;
1514  // set initial condition
1515  eq_data_->full_solution.add(loc_dofs[edge_dof], init_value/n_sides_of_edge);
1516  }
1517  }
1518 
1519  eq_data_->full_solution.ghost_to_local_begin();
1520  eq_data_->full_solution.ghost_to_local_end();
1521  eq_data_->full_solution.local_to_ghost_begin();
1522  eq_data_->full_solution.local_to_ghost_end();
1523 }
1524 
1525 
1527 {
1528  START_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
1529 
1530  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1531  unsigned int dim = dh_cell.dim();
1532  assembler[dim-1]->assemble_reconstruct(dh_cell);
1533  }
1534 
1535  eq_data_->full_solution.local_to_ghost_begin();
1536  eq_data_->full_solution.local_to_ghost_end();
1537 }
1538 
1539 
1541  // save diagonal of steady matrix
1542  MatGetDiagonal(*( schur0->get_matrix() ), steady_diagonal);
1543  // save RHS
1544  VecCopy(*( schur0->get_rhs()), steady_rhs);
1545 
1546 
1547  PetscScalar *local_diagonal;
1548  VecGetArray(new_diagonal,& local_diagonal);
1549 
1550  DebugOut().fmt("Setup with dt: {}\n", time_->dt());
1551 
1552  balance_->start_mass_assembly(eq_data_->water_balance_idx);
1553 
1554  std::vector<LongIdx> dofs;
1555  dofs.reserve(eq_data_->dh_->max_elem_dofs());
1556 
1557  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1558  ElementAccessor<3> ele = dh_cell.elm();
1559  dofs.resize(dh_cell.n_dofs());
1560  dh_cell.get_dof_indices(dofs);
1561 
1562  // TODO: replace with DHCell getter when available for FESystem component
1563  const uint p_ele_dof = dh_cell.n_dofs() / 2;
1564  // set new diagonal
1565  double diagonal_coeff = eq_fields_->cross_section.value(ele.centre(), ele)
1566  * ( eq_fields_->storativity.value(ele.centre(), ele)
1567  +eq_fields_->extra_storativity.value(ele.centre(), ele))
1568  * ele.measure();
1569  local_diagonal[dh_cell.get_loc_dof_indices()[p_ele_dof]]= - diagonal_coeff / time_->dt();
1570 
1571  balance_->add_mass_values(eq_data_->water_balance_idx, dh_cell,
1572  {dh_cell.get_loc_dof_indices()[p_ele_dof]},
1573  {diagonal_coeff}, 0.0);
1574  }
1575  VecRestoreArray(new_diagonal,& local_diagonal);
1576  MatDiagonalSet(*( schur0->get_matrix() ), new_diagonal, ADD_VALUES);
1577 
1578  schur0->set_matrix_changed();
1579 
1580  balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1581 }
1582 
1584  START_TIMER("modify system");
1585  if (time_->is_changed_dt() && time_->step().index()>0) {
1586  double scale_factor=time_->step(-2).length()/time_->step().length();
1587  if (scale_factor != 1.0) {
1588  // if time step has changed and setup_time_term not called
1589  MatDiagonalSet(*( schur0->get_matrix() ),steady_diagonal, INSERT_VALUES);
1590 
1591  VecScale(new_diagonal, time_->last_dt()/time_->dt());
1592  MatDiagonalSet(*( schur0->get_matrix() ),new_diagonal, ADD_VALUES);
1594  }
1595  }
1596 
1597  // modify RHS - add previous solution
1598  VecPointwiseMult(*( schur0->get_rhs()), new_diagonal, previous_solution);
1599 // VecPointwiseMult(*( schur0->get_rhs()), new_diagonal, schur0->get_solution());
1600  VecAXPY(*( schur0->get_rhs()), 1.0, steady_rhs);
1602 
1603  //VecSwap(previous_solution, schur0->get_solution());
1604 }
1605 
1606 
1607 /// Helper method fills range (min and max) of given component
1608 void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component) {
1609  if (component==0) {
1610  min = 0;
1611  max = n_dofs/2;
1612  } else if (component==1) {
1613  min = n_dofs/2;
1614  max = (n_dofs+1)/2;
1615  } else {
1616  min = (n_dofs+1)/2;
1617  max = n_dofs;
1618  }
1619 }
1620 
1621 
1623  ASSERT_LT(component, 3).error("Invalid component!");
1624  unsigned int i, n_dofs, min, max;
1625  std::vector<int> dof_vec;
1626  std::vector<LongIdx> dof_indices(eq_data_->dh_->max_elem_dofs());
1627  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1628  n_dofs = dh_cell.get_dof_indices(dof_indices);
1629  dofs_range(n_dofs, min, max, component);
1630  for (i=min; i<max; ++i) dof_vec.push_back(dof_indices[i]);
1631  }
1632  return dof_vec;
1633 }
1634 
1635 
1636 //-----------------------------------------------------------------------------
1637 // vim: set cindent:
LimitSide::right
@ right
Element::n_neighs_vb
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:65
LocDofVec
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
DarcyMH::balance_
std::shared_ptr< Balance > balance_
Definition: darcy_flow_mh.hh:365
Side::edge
Edge edge() const
Returns pointer to the edge connected to the side.
Definition: accessors_impl.hh:227
EquationBase::mesh_
Mesh * mesh_
Definition: equation.hh:219
linsys_BDDC.hh
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
Input::Type::Bool
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
result_zeros
@ result_zeros
Definition: field_algo_base.hh:74
LinSys::get_solution
const Vec & get_solution()
Definition: linsys.hh:281
MPI_MIN
#define MPI_MIN
Definition: mpi.h:198
DarcyMH::previous_solution
Vec previous_solution
Definition: darcy_flow_mh.hh:389
DarcyMH::solution_precision
virtual double solution_precision() const
Definition: darcy_flow_mh.cc:842
DarcyMH::init_eq_data
void init_eq_data()
Definition: darcy_flow_mh.cc:426
ASSERT_GT
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
LinSys::SolveInfo::n_iterations
int n_iterations
Definition: linsys.hh:109
Armor::vec
ArmaVec< double, N > vec
Definition: armor.hh:885
DarcyMH::solve_nonlinear
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
Definition: darcy_flow_mh.cc:703
DarcyFlowMHOutput::output
void output()
Calculate values for output.
Definition: darcy_flow_mh_output.cc:246
TimeGovernor::dt
double dt() const
Definition: time_governor.hh:565
DarcyMH::use_steady_assembly_
bool use_steady_assembly_
Definition: darcy_flow_mh.hh:374
UnitSI::dimensionless
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
SchurComplement::set_complement
void set_complement(LinSys_PETSC *ls)
Set complement LinSys object.
Definition: schur.cc:287
factory.hh
field_constant.hh
LinSys
Definition: la_linsys_new.hh:169
vector_mpi.hh
time_governor.hh
Basic time management class.
field_add_potential.hh
DarcyMH::size
int size
Definition: darcy_flow_mh.hh:369
DarcyMH::steady_rhs
Vec steady_rhs
Definition: darcy_flow_mh.hh:387
field_algo_base.hh
DarcyMH::solve_time_step
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
Definition: darcy_flow_mh.cc:642
MeshBase::region_db
const RegionDB & region_db() const
Definition: mesh.h:175
DarcyMH::update_solution
void update_solution() override
Definition: darcy_flow_mh.cc:627
DarcyMH::EqFields::EqFields
EqFields()
Creation of all fields.
Definition: darcy_flow_mh.cc:189
DarcyMH::EqFields::river
@ river
Definition: darcy_flow_mh.hh:154
LinSys::SolveInfo
Definition: linsys.hh:104
ASSERT
#define ASSERT(expr)
Definition: asserts.hh:351
DarcyMH::type_field_descriptor
static const Input::Type::Record & type_field_descriptor()
Definition: darcy_flow_mh.cc:119
LinSys::set_symmetric
void set_symmetric(bool flag=true)
Definition: linsys.hh:560
DarcyMH::postprocess
virtual void postprocess()
postprocess velocity field (add sources)
Definition: darcy_flow_mh.cc:792
Mesh::n_sides
unsigned int n_sides() const
Definition: mesh.cc:308
DarcyFlowInterface::MortarP1
@ MortarP1
Definition: darcy_flow_interface.hh:33
EquationBase::time_
TimeGovernor * time_
Definition: equation.hh:220
DarcyMH::setup_time_term
virtual void setup_time_term()
Definition: darcy_flow_mh.cc:1540
distribution.hh
Support classes for parallel programing.
LinSys_BDDC::get_input_type
static const Input::Type::Record & get_input_type()
Definition: linsys_BDDC.cc:39
bc_mesh.hh
Edge::n_sides
unsigned int n_sides() const
Returns number of sides aligned with the edge.
Definition: accessors.hh:329
Input::Type::Integer
Class for declaration of the integral input data.
Definition: type_base.hh:483
Input::Record::val
const Ret val(const string &key) const
Definition: accessors_impl.hh:31
DHCellAccessor::get_dof_indices
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
Definition: dh_cell_accessor.hh:80
TimeGovernor::set_upper_constraint
int set_upper_constraint(double upper, std::string message)
Sets upper constraint for the next time step estimating.
Definition: time_governor.cc:555
EquationBase::time
TimeGovernor & time()
Definition: equation.hh:148
Input::Type::Selection::close
const Selection & close() const
Close the Selection, no more values can be added.
Definition: type_selection.cc:65
fmt::fprintf
FMT_FUNC int fprintf(std::ostream &os, CStringRef format, ArgList args)
Definition: ostream.cc:56
FilePath
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
linsys.hh
Wrappers for linear systems based on MPIAIJ and MATIS format.
Input::Type::Double
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
LinSys::rhs_zero_entries
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:272
chkerr
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142
DarcyMH::create_linear_system
void create_linear_system(Input::AbstractRecord rec)
Definition: darcy_flow_mh.cc:1028
FLOW123D_FORCE_LINK_IN_CHILD
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
LinSys_PETSC
Definition: linsys_PETSC.hh:43
MeshBase::n_edges
unsigned int n_edges() const
Definition: mesh.h:114
LinSys::mat_zero_entries
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:263
THROW
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
LinSys::SolveInfo::converged_reason
int converged_reason
Definition: linsys.hh:108
std::vector< double >
ElementAccessor< 3 >
DarcyMH::output_data
virtual void output_data() override
Write computed fields.
Definition: darcy_flow_mh.cc:826
DarcyMH::~DarcyMH
virtual ~DarcyMH() override
Definition: darcy_flow_mh.cc:1447
system.hh
arma::vec3
Definition: doxy_dummy_defs.hh:17
LinSys::set_rhs_changed
void set_rhs_changed()
Definition: linsys.hh:217
fn_mh_velocity
Definition: assembly_models.hh:29
SchurComplement::make_complement_distribution
Distribution * make_complement_distribution()
get distribution of complement object if complement is defined
Definition: schur.cc:294
assembly_models.hh
Functors of FieldModels used in Darcy flow module.
DHCellAccessor::dim
unsigned int dim() const
Return dimension of element appropriate to cell.
Definition: dh_cell_accessor.hh:101
DarcyMH::get_mh_mortar_selection
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
Definition: darcy_flow_mh.cc:81
LinSys::solve
virtual SolveInfo solve()=0
field_fe.hh
uint
unsigned int uint
Definition: mh_dofhandler.hh:101
Input::AbstractRecord::type
Input::Type::Record type() const
Definition: accessors.cc:273
linsys_PETSC.hh
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
TimeStep::use_fparser_
bool use_fparser_
Definition: time_governor.hh:235
DarcyMH::max_n_it_
unsigned int max_n_it_
Definition: darcy_flow_mh.hh:380
dofs_range
void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component)
Helper method fills range (min and max) of given component.
Definition: darcy_flow_mh.cc:1608
LinSys::set_negative_definite
void set_negative_definite(bool flag=true)
Definition: linsys.hh:587
MeshBase::elements_range
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
index_types.hh
LinSys::set_positive_definite
void set_positive_definite(bool flag=true)
Definition: linsys.hh:575
ASSERT_LT
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
DarcyMH::allocate_mh_matrix
void allocate_mh_matrix()
Definition: darcy_flow_mh.cc:878
LinSys::set_matrix_changed
void set_matrix_changed()
Definition: linsys.hh:211
fe_p.hh
Definitions of basic Lagrangean finite elements with polynomial shape functions.
DarcyMH::prepare_new_time_step
void prepare_new_time_step()
Definition: darcy_flow_mh.cc:787
DarcyMH::assembly_mh_matrix
void assembly_mh_matrix(MultidimAssembly &assembler)
Definition: darcy_flow_mh.cc:854
LinSys::start_allocation
virtual void start_allocation()
Definition: linsys.hh:332
DarcyMH::solved_time
virtual double solved_time() override
Definition: darcy_flow_mh.cc:403
DarcyFlowMHOutput::get_input_type
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
Definition: darcy_flow_mh_output.cc:62
Input::Type::Record::size
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
SchurComplement
Definition: schur.hh:64
TimeGovernor::is_end
bool is_end() const
Returns true if the actual time is greater than or equal to the end time.
Definition: time_governor.hh:595
SchurComplement::set_from_input
void set_from_input(const Input::Record in_rec) override
Definition: schur.cc:138
Input::Type::Default
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
DHCellSide
Side accessor allows to iterate over sides of DOF handler cell.
Definition: dh_cell_accessor.hh:176
Input::Type::Record::derive_from
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
DarcyMH::eq_data_
std::shared_ptr< EqData > eq_data_
Definition: darcy_flow_mh.hh:392
Distribution
Definition: distribution.hh:50
DarcyMH::output_object
DarcyFlowMHOutput * output_object
Definition: darcy_flow_mh.hh:367
Mesh::bc_mesh
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
Definition: mesh.h:567
AddPotentialFactory
Definition: field_add_potential.hh:49
LimitSide::left
@ left
LinSys::compute_residual
virtual double compute_residual()=0
DarcyFlowInterface::NoMortar
@ NoMortar
Definition: darcy_flow_interface.hh:31
accessors.hh
LinSys::set_tolerances
virtual void set_tolerances(double r_tol, double a_tol, unsigned int max_it)=0
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
LinSys::rhs_set_value
void rhs_set_value(int row, double val)
Definition: linsys.hh:380
FieldFlag::equation_result
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
sys_profiler.hh
DarcyMH::EqFields::get_bc_type_selection
static const Input::Type::Selection & get_bc_type_selection()
Return a Selection corresponding to enum BC_Type.
Definition: darcy_flow_mh.cc:90
darcy_flow_mh_output.hh
Output class for darcy_flow_mh model.
ASSERT_PERMANENT
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
field_model.hh
DarcyMH::assembly_source_term
virtual void assembly_source_term()
Source term is implemented differently in LMH version.
Definition: darcy_flow_mh.cc:991
DarcyMH::DarcyFlowMHOutput
friend class DarcyFlowMHOutput
Definition: darcy_flow_mh.hh:394
LinSys::get_matrix
virtual const Mat * get_matrix()
Definition: linsys.hh:186
mpi.h
mixed_mesh_intersections.hh
schur.hh
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
DarcyMH::modify_system
void modify_system()
Definition: darcy_flow_mh.cc:1583
DarcyMH::EqFields
Definition: darcy_flow_mh.hh:143
DarcyMH::min_n_it_
unsigned int min_n_it_
Definition: darcy_flow_mh.hh:379
TimeGovernor
Basic time management functionality for unsteady (and steady) solvers (class Equation).
Definition: time_governor.hh:317
TimeGovernor::step
const TimeStep & step(int index=-1) const
Definition: time_governor.cc:756
DarcyMH::zero_time_term
virtual bool zero_time_term(bool time_global=false)
Definition: darcy_flow_mh.cc:694
DarcyMH::steady_diagonal
Vec steady_diagonal
Definition: darcy_flow_mh.hh:386
DHCellAccessor::elm
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
Definition: dh_cell_accessor.hh:71
field_values.hh
Input::AbstractRecord
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
Input::Type::Default::obligatory
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
SchurComplement::IsB
IS IsB
Definition: schur.hh:138
UnitSI
Class for representation SI units of Fields.
Definition: unit_si.hh:40
LinSys_BDDC
Definition: linsys_BDDC.hh:39
fn_mh_piezohead
Definition: assembly_models.hh:37
ASSERT_EQ
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
DarcyMH::print_matlab_matrix
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
Definition: darcy_flow_mh.cc:1197
EquationBase::input_record_
Input::Record input_record_
Definition: equation.hh:221
RegionDB::get_region_set
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
FieldCommon::field_descriptor_record_description
static const std::string field_descriptor_record_description(const string &record_name)
Definition: field_common.cc:73
DarcyMH::eq_fields_
std::shared_ptr< EqFields > eq_fields_
Definition: darcy_flow_mh.hh:391
Input::Type::Record::declare_key
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
MPI_INT
#define MPI_INT
Definition: mpi.h:160
DarcyFlowMHOutput::get_input_type_specific
static const Input::Type::Instance & get_input_type_specific()
Definition: darcy_flow_mh_output.cc:69
DarcyMH::schur0
LinSys * schur0
Definition: darcy_flow_mh.hh:384
mesh.h
DHCellAccessor::n_dofs
unsigned int n_dofs() const
Return number of dofs on given cell.
Definition: dh_cell_accessor.hh:420
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:65
std::map
Definition: doxy_dummy_defs.hh:11
DarcyMH::reconstruct_solution_from_schur
void reconstruct_solution_from_schur(MultidimAssembly &assembler)
Definition: darcy_flow_mh.cc:1526
DarcyFlowInterface::MortarMethod
MortarMethod
Type of experimental Mortar-like method for non-compatible 1d-2d interaction.
Definition: darcy_flow_interface.hh:30
DarcyFlowInterface::get_input_type
static Input::Type::Abstract & get_input_type()
Definition: darcy_flow_interface.hh:23
LinSys::start_add_assembly
virtual void start_add_assembly()
Definition: linsys.hh:340
DarcyMH::EqFields::none
@ none
Definition: darcy_flow_mh.hh:150
TimeGovernor::estimate_dt
double estimate_dt() const
Estimate choice of next time step according to actual setting of constraints.
Definition: time_governor.cc:631
DarcyMH::n_schur_compls
int n_schur_compls
Definition: darcy_flow_mh.hh:370
TimeGovernor::view
void view(const char *name="") const
Definition: time_governor.cc:772
FieldFlag::input_copy
static constexpr Mask input_copy
Definition: field_flag.hh:44
Input::Type::Record::close
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
LinSys_BDDC::set_from_input
void set_from_input(const Input::Record in_rec) override
Definition: linsys_BDDC.cc:279
SchurComplement::loc_size_B
int loc_size_B
Definition: schur.hh:137
Input::Type
Definition: balance.hh:41
MixedMeshIntersections::element_intersections_
std::vector< std::vector< ILpair > > element_intersections_
Definition: mixed_mesh_intersections.hh:83
partitioning.hh
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:182
DarcyMH::new_diagonal
Vec new_diagonal
Definition: darcy_flow_mh.hh:388
DarcyMH::get_component_indices_vec
std::vector< int > get_component_indices_vec(unsigned int component) const
Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
Definition: darcy_flow_mh.cc:1622
LongIdx
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
LinSys_PETSC::set_from_input
void set_from_input(const Input::Record in_rec) override
Definition: linsys_PETSC.cc:470
Mesh
Definition: mesh.h:362
DarcyMH::EqFields::total_flux
@ total_flux
Definition: darcy_flow_mh.hh:152
OutputTime::get_input_type
static const Input::Type::Record & get_input_type()
The specification of output stream.
Definition: output_time.cc:38
Element::n_sides
unsigned int n_sides() const
Definition: elements.h:131
DarcyMH::nonlinear_iteration_
unsigned int nonlinear_iteration_
Definition: darcy_flow_mh.hh:381
Input::Type::Record::copy_keys
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
FieldAlgorithmBase
Definition: field_algo_base.hh:112
Input::Type::Array
Class for declaration of inputs sequences.
Definition: type_base.hh:339
DarcyMH::tolerance_
double tolerance_
Definition: darcy_flow_mh.hh:378
DarcyMH::set_mesh_data_for_bddc
void set_mesh_data_for_bddc(LinSys_BDDC *bddc_ls)
Definition: darcy_flow_mh.cc:1270
DarcyMH::get_input_type
static const Input::Type::Record & get_input_type()
Definition: darcy_flow_mh.cc:134
Model
Definition: field_model.hh:338
DHCellAccessor
Cell accessor allow iterate over DOF handler cells.
Definition: dh_cell_accessor.hh:43
LinSys::finish_assembly
virtual void finish_assembly()=0
LinSys::get_rhs
virtual const Vec * get_rhs()
Definition: linsys.hh:202
MPI_Allreduce
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
Definition: mpi.h:612
Input::Array
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
DarcyMH::assembly_linear_system
virtual void assembly_linear_system()
Definition: darcy_flow_mh.cc:1141
DarcyMH::DarcyMH
DarcyMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
Definition: darcy_flow_mh.cc:351
WarningOut
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
ElementAccessor::region
Region region() const
Definition: accessors.hh:198
EquationBase::record_template
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
Definition: equation.cc:35
LinSys::get_solution_precision
virtual double get_solution_precision()=0
MixedPtr
Definition: mixed.hh:247
SchurComplement
SchurComplement SchurComplement
Definition: linsys.hh:91
TimeGovernor::is_default
bool is_default()
Definition: time_governor.hh:388
darcy_flow_mh.hh
mixed-hybrid model of linear Darcy flow, possibly unsteady.
LinSys::set_solution
void set_solution(Vec sol_vec)
Definition: linsys.hh:289
DarcyMH::initialize_specific
virtual void initialize_specific()
Definition: darcy_flow_mh.cc:578
DarcyMH::EqFields::seepage
@ seepage
Definition: darcy_flow_mh.hh:153
std::vector::data
T data
Definition: doxy_dummy_defs.hh:7
LinSys::get_input_type
static Input::Type::Abstract & get_input_type()
Definition: linsys.cc:29
ElementAccessor::idx
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
FilePath::output_file
@ output_file
Definition: file_path.hh:69
EquationBase::eq_fieldset_
FieldSet * eq_fieldset_
Definition: equation.hh:228
TimeGovernor::is_changed_dt
bool is_changed_dt() const
Definition: time_governor.hh:529
DarcyFlowInterface::MortarP0
@ MortarP0
Definition: darcy_flow_interface.hh:32
balance.hh
local_to_global_map.hh
DarcyMH::data_changed_
bool data_changed_
Definition: darcy_flow_mh.hh:375
MeshBase::n_elements
unsigned int n_elements() const
Definition: mesh.h:111
TimeStep::index
unsigned int index() const
Definition: time_governor.hh:159
IntersectionLocalBase::bulk_ele_idx
unsigned int bulk_ele_idx() const
Returns index of bulk element.
Definition: intersection_local.hh:77
DarcyMH::EqFields::dirichlet
@ dirichlet
Definition: darcy_flow_mh.hh:151
Mesh::mixed_intersections
MixedMeshIntersections & mixed_intersections()
Definition: mesh.cc:849
DebugOut
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
Input::Type::Selection::add_value
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.
Definition: type_selection.cc:50
ElementAccessor::side
SideIter side(const unsigned int loc_index)
Definition: accessors_impl.hh:131
TimeGovernor::next_time
void next_time()
Proceed to the next time according to current estimated time step.
Definition: time_governor.cc:670
TimeStep::length
double length() const
Definition: time_governor.hh:160
Input::Type::Default::optional
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
DarcyMH::read_initial_condition
virtual void read_initial_condition()
Definition: darcy_flow_mh.cc:1502
RegionIdx::bulk_idx
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Definition: region.hh:90
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:115
FESystem
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:101
DarcyMH::registrar
static const int registrar
Registrar of class to factory.
Definition: darcy_flow_mh.hh:400
DarcyMH::initialize
void initialize() override
Definition: darcy_flow_mh.cc:472
LinSys_PETSC::get_input_type
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
DarcyMH::zero_time_step
void zero_time_step() override
Definition: darcy_flow_mh.cc:581
field.hh
DarcyFlowInterface
Definition: darcy_flow_interface.hh:15
Balance::get_input_type
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:53
IntersectionResult::none
@ none
assembly_mh_old.hh
ElementAccessor::measure
double measure() const
Computes the measure of the element.
Definition: accessors_impl.hh:59
END_TIMER
#define END_TIMER(tag)
Ends a timer with specified tag.
Definition: sys_profiler.hh:149
DarcyMH::eq_fields
EqFields & eq_fields()
Definition: darcy_flow_mh.hh:268
TimeGovernor::last_dt
double last_dt() const
Definition: time_governor.hh:548
FieldValue
Definition: field_values.hh:645
range_wrapper.hh
Implementation of range helper class.
ElementAccessor::centre
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
Definition: accessors_impl.hh:81
IntersectionLocalBase
Common base for intersection object.
Definition: intersection_local.hh:48
Distribution::begin
unsigned int begin(int proc) const
get starting local index
Definition: distribution.hh:109
DarcyMH::EqData::EqData
EqData()
Constructor.
Definition: darcy_flow_mh.cc:331
MessageOut
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
LinSys::mat_set_values
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
intersection_local.hh
Classes with algorithms for computation of intersections of meshes.
TimeGovernor::t
double t() const
Definition: time_governor.hh:542